{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named 'statsmodels'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-f0600a0351ef>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinear_model\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mLinearRegression\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mstatsmodels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mapi\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0msm\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mstatsmodels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformula\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mapi\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0msmf\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mImportError\u001b[0m: No module named 'statsmodels'"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.arange(1, 25, 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# linear relationship with contant variance of residual\n",
    "# x1 和 y具有恒定的线性关系, 具有稳定的残差方差\n",
    "x1 = y.copy() + np.random.randn(96) # 生成96个正太分布的随机数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# non contant variance with residuals \n",
    "# x2和y2具有不稳定的残差方差\n",
    "x2 = y.copy()\n",
    "y2 = x2 + np.concatenate((np.random.randn(20)*0.5,\n",
    "                                np.random.randn(20)*1, \n",
    "                                np.random.randn(20)*4, \n",
    "                               np.random.randn(20)*6, \n",
    "                               np.random.randn(16)*8), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGECAYAAAAr/IHVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYXVWZ7/vfm1BAgUiBBDqUxKTVjrs1baLZ6NlpFYMSxe5NDN5oG8G2jd2nfVpoTFt6+gjtLbWNgJ6j2zaKG7yDApEmtoFjpNmybTQxQUDM9sK1iFyE4lpKEd7zx5wrWbVqznWd9/n9PE89qVprrrXGWlWZY75jvOMd5u4CAAAAAAxuTt4NAAAAAICqIMACAAAAgIQQYAEAAABAQgiwAAAAACAhBFgAAAAAkBACLAAAAABICAEWkDAz+4CZfSHvdnTLzFaY2S/M7FEzW53i6ywIX2NuzP3nmNlXEnotN7PnJPFcAICAmf2bmZ2Wdzu6ZWZ/a2b3hH3PM1J8nbea2VVt7r/GzP46gdc5zszuGvR5kD4CLPTNzP7CzLaFJ67d4Yn3TxN43gvN7CNJtDF8vkRObN1y94+5e2avl4APSfq0uz/N3Tel9SLufkf4GnvSeg0A1WVmt4UXywc33fbXZnZNjs3qSfgeXpXg82U6mOTur3X3i7J6vUGY2ZCk8ySdEPY9v03rtdz9q+5+QlrPj/IhwEJfzOwfJH1S0sckHSVpgaT/LumkPNuVNzPbL+829OFZkm7u5sCSvj8A1bGfpPfk3Yi6sUDZrhmPknSg6N+Qg7L9Z0EBmNmhCmY9/s7dL3P3x9x92t3/1d3XhcccYGafNLO7w69PmtkB4X3HmdldZnaWmd0bzn69PbxvraS3SvrHcGbsX8Pbx8zsV2b2iJn9zMxe39Se083sB2b2CTN70MxuNbPXhvd9VNLLJH06fL5PR7yf75rZu1tuu8HM1oTff8rM7jSzh81su5m9rOm4c8zsW2b2FTN7WNLpraluZvZNM/uNmT1kZtea2fOb7rvQzD5jZpvD93a9mT276f7nm9nVZvZAOHL7gfD2OU2fyW/N7BIzO7zN7+ydZvbL8HmuMLOjw9t/JekPJf1r+PkcEPHY28zsfWb2U0mPmdl+Zna0mV1qZveFn/ffNx1/bDiz+XDY5vPC2xeGo637hT8vMrN/D9/31ZKOaHqOWWkQzSO/4Wv80Mwmw7+fT5vZ/jHv/cTwb+YRM5sws/fGfU4ACm+DpPea2UjUnWb2X8zsx+H59sdm9l+a7rvGzD5sZteF54OrzOyIqOcJjz/JzHaG57JfmdlrwtuPDs+jD4Tn1Xc2Peac8Hz8pfA1bjaz5eF9X1YwGNk43/5jeHtffYSZXRsedkP4fG9uaf8B4TnyBU23zTOzKTM70swOM7Mrw/P4g+H3z2z5vD5qZtdJelzSH1pTRoiZPdvMtoZ90P1m9tXm30t4zn6vmf00fG8Xm9mBXXy+h5rZBeG5fcLMPmLxqeWR1xpm9keSdoWHTZrZ1ojHNvqkd5jZHZK2hre/1Mz+V/jZ3WBmxzU95nQz+3X4u7jVzN7adPsPmo57tZn9PHzfn5ZkLX8jX4loR6NvfLuZ3RK+xq/N7F1R7z089n3hZ/SIme0ys+PjjkXG3J0vvnr6kvQaSU9K2q/NMR+S9B+SjpQ0T9L/kvTh8L7jwsd/SNKQpBMVnLwPC++/UNJHWp7vjZKOVjAo8GZJj0maH953uqRpSe+UNFfS30q6W5KF918j6a/btPVtkq5r+vmPJU1KOiD8+S8lPUPByOlZkn4j6cDwvnPC114dtm04vO0rTc/3V5IOkXSAglm/nU33XSjpAUnHhs//VUnfCO87RNLu8DUPDH9+SXjfGeHn+8zweT8n6esx72+lpPslvSg89v+VdG3T/bdJelWbz+c2STslHRO+vzmStkv6oKT9FQRov5a0Kjz+h5JODb9/mqSXht8vlOSNv5vwuPPCNr1c0iONzy38G7kroh2vCr9/saSXhp/ZQkm3SDqj6ViX9Jzw+92SXhZ+f5ikF+X9f4gvvvjq/atxDpB0mcI+QtJfS7om/P5wSQ9KOjU8N5wS/vyM8P5rJP1K0h+F57JrJI3HvNaxkh6S9OrwnDcq6Xnhff+uIGPjQElLJd0n6fjwvnMk/U5BvzZX0npJ/9H6Hlpeq68+Irx/77ku5n18UdJHm37+O0nfDb9/hqSTJR0Uvv43JW1qOvYaSXdIen742kNq6k8lPSf8fA5Q0M9fK+mTLe/1Rwr67sPD8/TfdPH5blLQpx2s4BriR5LeFfP+2l1rLFRTnxPx2Mb9Xwpfazhsx2/D39+csH2/DZ/7YEkPS1ocPn6+pOeH358u6Qfh90eEx70h/MzOVHDN0/jcztHMa4QZ7ZT0OknPVhCUvULB9dGLwvuOU9g3Slos6U5JRzc9z7Pz/n/KV/h7zbsBfJXvS8EM0286HPMrSSc2/bxK0m3h98dJmmo+6Um6V/suxC9US4AV8fw7JZ0Ufn+6pF823XdQeLL6g/Dna9Q+wDpEQcD2rPDnj0r6YpvjH5T0wvD7c9QUrDTd9pWYx46EbTu06b1+oen+EyX9PPz+FEk7Yp7nFoUdevjzfAWB3qyORNIFkj7e9PPTwmMXhj/fps4B1l81/fwSSXe0HPN+Sf8j/P5aSf8s6YiWY/Z2IgpGcZ+UdHDT/V9TlwFWRBvPkHR508/NAdYdkt4l6el5/9/hiy+++v/SvgDrBQouzudpZoB1qqQftTzmh5JOD7+/RtI/Nd33fyoMNiJe63OSzo+4/RhJeyQd0nTbekkXht+fI+n/a7rvjyVNtb6HNu+x6z4i/LlTgPUqSb9u+vk6SW+LOXappAebfr5G0odajrlGMf2pgoHGHU0/3ybpL5t+/rikf+nw+R4l6feShptuO0XS92Nes921xkJ1F2D9YdNt75P05Zbjtkg6TUGANakgKB1uOeZ07Quw3qaZQbVJuktdBlgR7dwk6T3h98dpX4D1HAXXTq+SNJTF/0G+uv8iRRD9+K2kI6x9vvLRkm5v+vn28La9z+HuTzb9/LiCC/9IZva2MJVg0swmFXSwzakdv2l84+6Ph9/GPl8zd39E0mZJbwlveouCUcLGa58VTtc/FL72oS2vfWebds81s/Ew/eFhBR2O4tqumZ/DMQo6jyjPknR50+dxi4JO/6iIY2f8Ltz9UQW/w9G4dkdofo/PknR047XD1/9A02u/Q8EI8c8tSNH5s5g2PejujzXddnvEcZHM7I/CdJbfhJ/rxzTzM212soKLktstSEn8P7p9HQDF4+43SbpS0ljLXa39jsKfm891cefbVnHn36MlPRD2G92+xoFx/eWAfUQ3tkoaNrOXmNmzFARRl4evfZCZfc7Mbg9f+1pJIy3peO36tyPN7BthitrDkr6i2efhXvu3ZymY9dnd1L98TsEMVZRO1xrdaO3f3tjSv/2pgoyZxxRk0PxN2L7NZva8mDbtfU4PoqHYz7GVmb3WzP7DghTUSQX916z+zd1/qWBw8RxJ94a/i17fO1JCgIV+/FBBCkS7kt53KzhRNSwIb+uGN/8Qdgqfl/RuBakeI5JuUlNOcy/PF+Prkk4JL76HJX0/fO2XKRjRepOCFMYRBSOnza/d7vn/QkHhj1cpCMwWhrd30/Y7FaQJxN33Wncfafo60N0nIo6d8buwoALXMyRFHRun+T3eKenWltc+xN1PlCR3/4W7n6KgQ/xvkr5lTVW/QrslHdZy+4Km7x9TMBPZaPNcBaPVDZ+V9HNJz3X3pysI8CI/U3f/sbufFLZnk6RLun7XAIrqbAVp4c2BTWu/IwXnlV7OdQ1x59+7JR1uZof0+Rqt/cUgfUTnF3N/SsE575Twta5sCg7PUpBm9pLwPPryiNdu17+tD+//k/Dxf9lDu+M+3zsVzGAd0dS/PN3dnx9xrDTYtUZDa//25Zb+7WB3H5ckd9/i7q9WkDXycwXXJq12KwggJQUFQpp/Vkv/JukPmo49QNKlkj4h6ajwmuM7iu/fvubuf6rgM3AFfS4KgAALPXP3hxSsv/mMma0OR8GGwlGXj4eHfV3SP1mwoPaI8Phu9zi6R8G6noaDFZw47pOCBaAKZrC61fp8Ub6j4AT1IUkXh52SFKQPPhm+9n5m9kFJT+/htQ9R0Fn8VsEJ9WM9PPZKSX9gZmeEi3YPMbOXhPf9i6SPhsFnY+FyXAXHr0l6u5ktDU/eH5N0vbvf1kNbmv1I0sPh4trhcAT2BWb2n8O2/KWZzQs/w8nwMTNKs7v77ZK2SfpnM9vfgvL+f950yP9WMOr7OgtK7f6Tgjz/hkMU5Lg/Go4g/m1UQ8PnfquZHeru0+FjKBMPlFw4en+xpL9vuvk7kv7Igi1E9rOg6MMfKziX9uoCBefN4y0oKjRqZs9z9zsVrPNZb2YHmtmfKJi1/2rbZ9untT8apI+Ier4oX1Mw8/LW8Pvm155SUATicAVBay8OkfRo+PhRSet6eGzc57tb0lWSzjWzp4f3PdvMXhHzPINca0T5iqQ/N7NVYd92oAVFl55pZkeZ2X8NBwZ/r+C9R/UnmyU938zWhDOXf6+mIErBEoeXW7A35KEKUuwb9lfQ190n6UkLCnZFln83s8VmtjLs13+n4HdJ/1YQBFjoi7ufJ+kfFFz43qdg1OfdCmYIJOkjCi6gfyrpRkk/CW/rxgWS/jicnt/k7j+TdK6CmbN7JC1RkEferU9JeoMFVZL+n5j383sFC6dfpZkd0BZJ/6bggv92BSexrqf6FSyevV3B6ObPFCzG7Uo4yvhqBYHHbyT9QtIrw7s/JekKSVeZ2SPh874k5nm+J+n/VjAqtlvBqOFboo7tsl17wjYtlXSrggIaX1Aw+ioFRVBuNrNHw3a+xd1/F/FUfxG2+QEFHfuXml7jIQXrI76g4LN7TEEOe8N7w8c/omAE8eI2TT5V0m1hCsvfKBhlBVB+H1IwACdJ8mCfoz9TMDPzW0n/KOnP3P3+Xp/Y3X8k6e2SzleQtfDv2jdTcoqCmaa7FaTbne3uV3f51OsVBASTFlQ07buPCJ0j6aLw+d4U816uV3AOPVpBf9bwSQUZG/eHr/vdHl/7nxUUT3pIQVBxWbcP7PD5vk1BoPEzBWuev6VgxijKINcaUe26U8GM4ge079pmnYLr5TkK/rbuVtBvvUJBP9X6HPcrKMw1ruDv8LlqumYJ/1YuDtu8XU0DAGG///cKZh0fVNDPXRHT3APC17hfwTXCkWG7UQCNKmsAAAAAgAExgwUAAAAACSHAAgAAAICEEGABAAAAQEIIsAAAAAAgIe02ii2MI444whcuXJh3MwAAKdq+ffv97j6v85HFRF8FANXWbT9VigBr4cKF2rZtW97NAACkyMxuz7sNg6CvAoBq67afIkUQAAAAABJCgAUAAAAACSHAAgAAAICEEGABAAAAQEIIsAAAAAAgIQRYAAAAAJAQAiwAAAAASEhqAZaZHWNm3zezW8zsZjN7T3j7OWY2YWY7w68T02oDAAAAAGQpzY2Gn5R0lrv/xMwOkbTdzK4O7zvf3T+R4msDAAAAQOZSC7Dcfbek3eH3j5jZLZJG03o9AAAAAMhbmjNYe5nZQknLJF0vaYWkd5vZ2yRtUzDL9WDEY9ZKWitJCxYsyKKZAFBKm3ZMaMOWXbp7ckpHjwxr3arFWr2M8SwAQLFVtf9KvciFmT1N0qWSznD3hyV9VtKzJS1VMMN1btTj3H2juy939+Xz5s1Lu5kAUEqbdkzo/ZfdqInJKbmkickpvf+yG7Vpx0TeTQMAIFaV+69UAywzG1IQXH3V3S+TJHe/x933uPtTkj4v6dg02wAAVbZhyy5NTe+ZcdvU9B5t2LIrpxYBANBZlfuv1FIEzcwkXSDpFnc/r+n2+eH6LEl6vaSb0moDAFTd3ZNTPd2epKqmdgAA0pdn/5W2NNdgrZB0qqQbzWxneNsHJJ1iZksluaTbJL0rxTYAQKUdPTKsiYjO6OiR4VRft5Ha0Rh9bKR2SCpVkGVmB0q6VtIBCvrEb7n72WZ2oaRXSHooPPR0d98Z/SwAgF7l1X9lIbUUQXf/gbubu/+Juy8Nv77j7qe6+5Lw9v/aNJsFAOjRulWLNTw0d8Ztw0NztW7V4lRft0KpHb+XtNLdX6hgbfBrzOyl4X3rmvovgisASFBe/VcWMqkiCABIR2O2KOtUvaqkdri7S3o0/HEo/PL8WgQA9ZBX/5UFAiwAKKG81z9VKbXDzOZK2i7pOZI+4+7Xm9nfSvqomX1Q0vckjbn77/NsJwBUzeplo5UIqFqlXqYdAJCsIpS2rVJqR1jZdqmkZ0o61sxeIOn9kp4n6T9LOlzS+6Iea2ZrzWybmW277777MmszAKC4CLAAoGSKsP5p9bJRrV+zRKMjwzJJoyPDWr9mSalHIt19UtI1kl7j7rs98HtJ/0MxW4qwZyMAoBUpggBQQO1SAIuy/qkKqR1mNk/StLtPmtmwpFdJ+m+NLUXCLUdWiy1FAABdIsACgILpVAK9SuufCmC+pIvCdVhzJF3i7lea2dYw+DJJOyX9TZ6NBACUBwEWABRMuxTA1ctGtW7V4hkBmFTe9U95c/efSloWcfvKHJoDAKgAAiwAKJhOKYBRpW1f+bx52rBll868eGelSt0CAFA2BFgAUDDdpAA2r3/qlFIIAACyQxVBACiYXkugF6GqIAAACBBgAUDB9FoCPS6lcGJyKtO9sQAAACmCAFBIvZRAj0splESqIABgYO22DsFszGABQMlFpRQ2kCoIABhEY53vxOSUXPvW+ZIhEY8ZLAAoucYo4hkX74y8v5sNiBmdBABE6bR1CGZjBgsAKmD1slGNxmw03GkDYkYnAQBxOm0dgtkIsACgInqtPthAFUIAQJy4QbpOg3d1RoAFABURVX3w5BePasOWXVo0tlkrxrdGzkoxOgkAiNPv4F2dsQYLACqknw2Iu9nYGABQT43+gnW63SPAAoCK6nZh8rpVi2cEYhKjkwCAfXrZOgQEWABQGr1W+us29Y/RSQAAkkOABQAl0G26X7NeUv8YnQQAIBkUuQCAEuin0h8LkwEAyB4zWABQAJ3S//qp9EfqHwAA2SPAAoCcdZP+12+lP1L/AADIFimCAJCzbtL/SPcDAKAcmMECgJx1k/5Huh8AAOVAgAUACem1jHpDt+l/pPsBAFB8pAgCQAIa66gmJqfk2reOatOOiY6PJf0PAIDqIMACgAT0U0a9YfWyUa1fs0SjI8MySaMjw1q/ZgmzVQAAlBABFgAkoJ8y6gAAoHoIsAAgAXHl0juVUZcGSy8EAADFQoAFAAkYZB3VIOmFAACgWAiwAKAPm3ZMaMX4Vi0a26wV41slqe91VKQXAgBQHZRpB4AeNVL6GrNOjZS+9WuW6LqxlT0/X7dl2gEAQPExgwUAPUo6pY8y7QAAVAcBFgD0KC51b2JySivGt/ZcnIIy7QAAVAcpggDQo7iUPmlfuqCkngKk1ctGCagAAKgAZrAAoEdRKX3NqAAIAEB9MYMFADE27ZjQhi27dPfklI4eGda6VYtnzDRt2LIrdiaLCoAAANQTM1gAEKHT5r+rl43qurGVGh1gg2EAAFA9BFgAEKHbSoFUAAQAAM1IEQSACN1u/tucLtiaShiXYojiMLMDJV0r6QAFfeK33P1sM1sk6RuSDpf0E0mnuvsT+bUUAFAWBFgAEKGXzX+jKgDGbUbcOB6F8XtJK939UTMbkvQDM/s3Sf8g6Xx3/4aZ/Yukd0j6bJ4NBQCUAymCABBh0NS/pDcjRjo88Gj441D45ZJWSvpWePtFklbn0DwAQAkxgwWg1rqpFNhPil+3KYbIn5nNlbRd0nMkfUbSryRNuvuT4SF3SYr8xZvZWklrJWnBggXpNxYAUHgEWABqq1Ma3yCb//aSYoh8ufseSUvNbETS5ZL+U9RhMY/dKGmjJC1fvjzyGABAvZAiCKC20kzjo7pg+bj7pKRrJL1U0oiZNQYhnynp7rzaBQAoFwIsALWVZhrf6mWjWr9miUZHhmWSRkeGtX7NEgpcFIyZzQtnrmRmw5JeJekWSd+X9IbwsNMkfTufFgIAyoYUQQC1lXQaX9R6ruvGVg7aTKRrvqSLwnVYcyRd4u5XmtnPJH3DzD4iaYekC/JsJACgPAiwANTWulWLZ6zBkvpP46Msezm5+08lLYu4/deSjs2+RQCAsiNFEEBtJZnGR1l2AAAgMYMFoObiKgXGlW+PQ1l2AAAgMYMFALM00v0mJqfk2pfut2nHROxj4tZtUZYdAIB6IcACgBb9pPtRlh0AAEgppgia2TGSviTpDyQ9JWmju3/KzA6XdLGkhZJuk/Qmd38wrXYAqI9e0/ri9JPu13idJF4fAIA4SfV1ZVb0zyDNNVhPSjrL3X9iZodI2m5mV0s6XdL33H3czMYkjUl6X4rtAFADSVbx67d8e9x6LgAAkkDF2nJ8BqmlCLr7bnf/Sfj9Iwo2bhyVdJKki8LDLpK0Oq02AKiPJKv4ke4HACgiKtaW4zPIpIqgmS1UsM/I9ZKOcvfdUhCEmdmRMY9ZK2mtJC1YsCCLZgIosSSr+PWa7lf0VAUAQDVQsbYcn0HqAZaZPU3SpZLOcPeHzayrx7n7RkkbJWn58uWeXgsBVEG/aX1xuk33K0OqAgCgGpLu68qoDJ9BqlUEzWxIQXD1VXe/LLz5HjObH94/X9K9abYBQD3kldZXhlQFAEA1kMJejs8gzSqCJukCSbe4+3lNd10h6TRJ4+G/306rDQDqo5e0viRT+sqQqgAAqAYq1pbjM0gzRXCFpFMl3WhmO8PbPqAgsLrEzN4h6Q5Jb0yxDQBqpJu0vqRT+sqQqgAAqA4q1hb/M0gtwHL3H0iKW3B1fFqvCwDttEvp6+dkvW7V4hkBm1S8VAUAALJC4aeMqggCQJbandyTTukrQ6oCAABZoPBTgAALQKV0OrmnkdJX9FQFAACykHSWSFmlWkUQALLWqapfGaoPAQBQRhR+ChBgAaiMTTsmImenpH0n99XLRrV+zRKNjgzLJI2ODGv9miW1GlkDACANcdkgdSv8RIoggEpopAbGaT65k9IHAEDyKPwUIMACUAlRqYENdTy5AwCQNQo/BQiwABRWL6Ve2+V3kwIIAEA2yBJhDRaAgmqk/E1MTsm1rxrgph0TkcfH5XePjgzX/kQPAACyQ4AFoJA6VQNsRXVAAABQBKQIAiikXku9kvcNAACKgAALQCHFbQg8ctBQ7GPI+wYAAHkjRRBAIa1btVhDc23W7Y/+7snYdVgAAKCcNu2Y0IrxrVo0tlkrxreWuq8nwAJQSKuXjerg/WdPsk8/5bHrsAAAQPn0Wtiq6EgRBFBYD01NR95+9+RUTyXcAQBAcbUrbFXGvp0ZLACFFVd6feSgoUqNdAEAUGe9FrYqOgIsAIUVV3rdXT2VcAcAAMUVN6Aad3vREWABSEQai1NXLxvV+jVLNDoyLFOwafD6NUvapg4CAIByqdpelqzBAjCwxuLUxqxSI2VP0sC501Gl1zds2RVZwr2sI10AABRJ1uucq7aXJQEWgIFlvTh13arFMwI6qdwjXQAAFEWag6btVGkvS1IEAQws68Wpq5eN6uQXj2quBftkzTXTyS+uzokZAIC8tBs0RXcIsAAMLOvFqZt2TOjS7RPa4y5J2uOuS7dPUEUQAIABVa2iXx4IsAAMLOvFqYyuISlmdoyZfd/MbjGzm83sPeHt55jZhJntDL9OzLutAJCFqlX0ywNrsAAMrJvFqUkumGV0DQl6UtJZ7v4TMztE0nYzuzq873x3/0SObQOAzLHOeXAEWAAS0W5xatILZo8eGaaKIBLh7rsl7Q6/f8TMbpHEYj4AtVXEin5ZVzUcFAEWgNQlXWWQ0TWkwcwWSlom6XpJKyS928zeJmmbglmuByMes1bSWklasGBBZm0FgDQVqaJfXlUNB8EaLACJa910OGq2Seo/pS9uA+KinmhRfGb2NEmXSjrD3R+W9FlJz5a0VMEM17lRj3P3je6+3N2Xz5s3L7P2AkBdlHHdNTNYABIVNdJkkjzi2EFS+oo0uoZyM7MhBcHVV939Mkly93ua7v+8pCtzah4A1FoZ110TYAFIVNRIk0uzgixS+lAEZmaSLpB0i7uf13T7/HB9liS9XtJNebQPAMosibVTZVx3TYAFIFFxI0rNwdVhBw3p7D9/PjNQKIIVkk6VdKOZ7Qxv+4CkU8xsqYI/3dskvSuf5gGou7IVeGhIau1UGdddE2ABSFTcSFOz300/lVFrgPbc/QcKJlhbfSfrtgBAqzIWeGhIqsBVEasadkKABSBRUSNNrQapIAgAQF0kXYU3S0munSrbumuqCAJIVGuFvzhFXpwKAEARlLHAQ0PcGqkir51KCgEWgMStXjaq68ZW6tbx12k0gRNsa9n3TTsmkmoqAACFVeYgZd2qxRoemjvjtqKvnUoKARaAVA16gm3kn09MTsm1L/+cIAsAUHVlDlLqvGcla7AApGrQxallzj8HAGAQZSzw0Kxsa6eSQoAFIHWtHURj9/VuTrplzj8HAGBQdQ1SyowUQQCpGyTNr8z55wAAoH4IsACkrl2aXydlzj8HAAD1Q4oggIF0s8P8IGl+Zc8/BwAA9UKABaBv3e4wf/TIsCYigqlu0/zIPwcAAGVBiiCAvnWb+keaHwAAqAtmsAD0rdvUv6TS/LpJRwQAAMgTARaAvvWS+jdoml+36YgAAAB5IkUQQN+yTP0bpBIhAABAVpjBAtC31ctGte32B/T16+/UHnfNNdPJL06nIAUbDgMAyoS09vpiBgtA3zbtmNCl2ye0x12StMddl26f6GoD4V6x4TAAoCwaae0Tk1Ny7UtrT6N/RPEQYAHoW5Zpe1QiBACUBWnt9UaKIIBZuk1ryDJtjw2HAQBFE9dfktZebwRYAGbopVrfoBsI94oNhwEARdGuv8y6f0SxkCIIYIZe0hpI2wMAlN2mHRNaMb5Vi8Y2a8X41q7XSbXrL+kf640ZLAAzRI24SdFpDaTtAQDKbJA9FtulAdI/1hsBFoC9Nu2YkEnyiPvi0hpI2wMAlFW7WahOfVunNED6x/qKerbKAAAgAElEQVQiRRDAXhu27IoMrkwirQEAUDmDFKMgDRBxUguwzOyLZnavmd3UdNs5ZjZhZjvDrxPTen0A8eLyzeM6FFfnVAkAAMpmkD0WVy8b1fo1SzQ6MiyTNDoyrPVrltBfItUUwQslfVrSl1puP9/dP5Hi6wJoo5+qR6NUPQIAVNC6VYtn9IlSb7NQpAEiSmozWO5+raQH0np+AP2h6hEAAAFmoZCGPIpcvNvM3iZpm6Sz3P3BqIPMbK2ktZK0YMGCDJsHVBtVjwAA2IdZKCQt6wDrs5I+rGBJx4clnSvpr6IOdPeNkjZK0vLly6PW3QPoQusu8yMHDenBx6dnHUfVIwAAZmvtRxl4RCeZBljufk/jezP7vKQrs3x9oG6i1lsNzTENzTVN79k3bkEaIAAAsw2yT1bZEEgmJ9MAy8zmu/vu8MfXS7qp3fEABhO13mr6KdfI8JAOPmA/TqIAgFhccA+2T1aZ1CmQzEJqAZaZfV3ScZKOMLO7JJ0t6TgzW6ogRfA2Se9K6/UBxK+3emhqWjvPPkHSvg70zIt31rYDBQDMxAV3YJB9ssqkLoFkVlILsNz9lIibL0jr9QDM1mmXeTpQAEAULrgDnfrRqqhLIJmV1Mq0A8hfp7Lr7TpQAEB9ccEdKPv2JZt2TGjF+FYtGtusFeNbtWnHRORxg2y4jNkIsIAK67S/RxIdaLcnb6CIzOwYM/u+md1iZjeb2XvC2w83s6vN7Bfhv4fl3VYgS2W/4E6qb2rtR0eGh3Tg0BydefHOwvd5jSyVickpufZlqUS1ueyBZNHksQ8WgAy1K7s+aOoDKYaogCcV7Mn4EzM7RNJ2M7ta0umSvufu42Y2JmlM0vtybCeQqXWrFs84v0vlueBOum9q9KNl6/N6SfNkH8xkEWABNTZoB0qOPsourGy7O/z+ETO7RdKopJMUFGqSpIskXSMCLNRImS+40+qbytbn9ZqlUtV9MPOohkmABVRAvyePQTtQcvRRJWa2UNIySddLOqqxrYi77zazI3NsGpCLsl5wp9U3la3Pq0uBjnbymnUkwAJKbtCTxyAdKCdvVIWZPU3SpZLOcPeHzazbx62VtFaSFixYkF4DAXQtrb6pbH1emdM8k5LXrCNFLoCSy7MSIItiUQVmNqQguPqqu18W3nyPmc0P758v6d6ox7r7Rndf7u7L582bl02DAbSVVt9Utj6vU6GrOshr1pEZLKAg+k3zyzNlocw5+oAkWTBVdYGkW9z9vKa7rpB0mqTx8N9v59A8AH1Iq28qY59X1jTPpOQ160iABRTAIGl+eacs1P3kjdJbIelUSTea2c7wtg8oCKwuMbN3SLpD0htzah+APqTVN9HnlUteaZIEWEABDJIjTI410D93/4GkuAVXx2fZFgDpyaOSHPKX16wjARZQAO3S/Dp1CmVMWQAAIEoagVDZ9q9CsvKYdSTAAgogLs1v5KChrjoFUhYAAGWXViBUtv2rUH5UEQQKIK4ykbtyqxAIAECW0qqKW7b9q1B+BFhAAcSVUn1oajry+E6dwqYdE1oxvlWLxjZrxfhWbdoxkUKrAQBITlqBUFzRp6LuX4XyI0UQKIioNL8NW3b1XCGQXHMAQBmlVRWXYlDIGjNYQIH1s6lhnhsPAwDQr7Q28i3DhrtknlQLM1hAgfVTIZBccwBAGaVZFbfIxaDIPKkeAiwgY72WoO21U8h742EAAPrV2uc1ZnaqvA0JVQ6rhxRBIEONUaqJySm59o1SJZkKkFaKBQAAWcqizywCMk+qJzbAMrPvmNnC7JoCVF8W66PKkGsOJIW+CshP2uuG6rKmuGxVDlkv1lm7FMELJV1lZhdJ+ri7R9eLBtC1rEapotIKe01NBEriQtFXAZnLYt1QXWZ2ylTlkPVi3YkNsNz9EjPbLOmDkraZ2ZclPdV0/3kZtA+olLzWR3FCRFXRVwH5yGLdUJ59ZpYDkmkW90ga68W606nIxbSkxyQdIOkQNXVaAHqX1ygVJ0RUHH0VkLEsZpfy6DPzGpAscpXDZnWZVRxUbIBlZq+RdJ6kKyS9yN0fz6xVQEUNOkrV76gaJ0RUFX0VkI8sZpfymNlhQLI9KhV3p90M1v8l6Y3ufnNWjQHqoN9RqkFG1TghosLoq4AcZDW7lPXMDgOS7ZVpvVieYqsIuvvL6LCA4hikmhKl21FV9FVAPqpasbZsFf2yVtXfe9LYaBhIQRoLZAcZVSvTAloAQDmUZd1QL5ih6WyQ33tdKhoTYAEJS2uB7KBpflXsCAEASBIDkumpU0VjAiwgYWktkGVUDQCAfdKaDWFAMh1ZFRApwiwZARaQkMZ/6KhZJmnwBbKMqgEAEKjTbEhVZFFApCh/FwRYQAJa/0NHSWKBLKNqAABQTr2MsqhoXJS/i9gqggC6F/UfuhmpfAAAJIdy6uWTRUXjovxdMIMFJKDdf9zRDql8RcgVBgCgTNjfsXyyWOpQlL8LAiwgAXH/oUdHhnXd2MrYxxUlVxgAgDKpeuGnqg6+pr3UoSh/F6QIAgnod9p7kM2DAQCoqypveNsYfJ2YnJJr3+Drph0TeTet8Iryd8EMFpCAfqe941ILJyantGhsc6VGrQAASFJVCz/VqZx5Gorwd0GABSSkn//QcamFkmaMWjWeHwCAJFX1IrvM6lTOvKpIEQRyFJVa2IqUQQBAGkhFK6a4ggxZlTPH4AiwgBytXjaqk188qrlmbY+j7CwAIGl1vsjetGNCK8a3atHYZq0Y35pIUJnUc9apnHlVkSII5GjTjgldun1Ce9zbHkfZWQBA0up6kZ1GelySz1mncuZVRYAF5KjTBsVStcrOAgCKo64X2WkUkUj6OetSzryqSBEEctRulLBqZWcBAMWSRSpaEaUxc1e22cCilDOvKmawgCZZV1Pqd4NiAAAGlUUqWtKS6KfTmLkr42xgEcqZVxUBFhDKo2QpU/QAgDyV6SI7qX46jb6X/hzNSBEEQnlUU2KKHsiXmX3RzO41s5uabjvHzCbMbGf4dWKebQQQSKqfTqPvzas/T6MaIgbHDBYQKlv+NIBEXCjp05K+1HL7+e7+ieybAyBOkv10GjN3Wc8GsllwcTGDBYSy2NivFZs8Avly92slPZB3OwB0lkc/XWRxM3pnXXIDM1o5I8ACQnlUU6rzJo9Awb3bzH4aphAeFneQma01s21mtu2+++7Lsn1A7dS16mGcuJm7Pe4M2uaMFEGUUhrV/vKopkRaIlBIn5X0YUke/nuupL+KOtDdN0raKEnLly9vv2M4kKKsq+DmoYxVD9MUV7mwWfNeXHX4GykKAiyUTpo5x1nnT5exrCtQde5+T+N7M/u8pCtzbA7QUdnW4gxyoV+mqodpi6pcGOXuyanS/Y2UHSmCKJ0qpdWR7gAUj5nNb/rx9ZJuijsWKIIy9YusPU5Oa+XCuWaRxx09Mlyqv5EqIMBC6ZQ5ra61nKokyrQDOTKzr0v6oaTFZnaXmb1D0sfN7EYz+6mkV0o6M9dGAh2UqV/kQj9Zq5eN6rqxlbp1/HU6900vjB20LdPfSBWQIojSKWtaXdz0/Po1S3Td2MqcWwfUk7ufEnHzBZk3BBhAXL/oklaMby3UWhsu9NPTbo3ahi27SnntVFbMYKF0yppWx6gdACANUf1iQ9FS8Ci1nq7mGa3rxlbuDbrKeu1UVqkFWGFp23vN7Kam2w43s6vN7Bfhv7Glb4E4ae+Wntau6IzaAQDS0NwvRinSYB4X+vlI+9oJM6WZInihpE9L+lLTbWOSvufu42Y2Fv78vhTbgIpKq4pQmlV2ypraCAAovka/uGhss6L2CyjKYB6l1vNDBcbspBZgufu1Zraw5eaTJB0Xfn+RpGtEgIUCaZfGN+hJKaqcKqN2AIAklWEwjwt9VF3Wa7COcvfdkhT+e2TcgWa21sy2mdm2++67L7MGop4aaYFxG/YlMfLH9DwAIG2k4AH5K2wVQXffKGmjJC1fvjxqthtIRGtaYJSkRv4YtQMApIkUPCB/WQdY95jZfHffHW7keG/Grw/MEpUW2IyRPwBAmTCYB+Qr6wDrCkmnSRoP//12xq+PGtm0Y6KrEbx26X+jjPwBAAAUWrfXfFlJLcAys68rKGhxhJndJelsBYHVJWb2Dkl3SHpjWq+PeuulGmDcguDRkWE2AAYAIGVFuzhGuaRZAbpfqRW5cPdT3H2+uw+5+zPd/QJ3/627H+/uzw3/fSCt10e99bKpLwuCAQDIR+PieGJySq7ibYyclrT23KyjXq75spJ1FUEgE71s6kt1PwAA8lHEi+O01TWoTEsv13xZIcBCJcVV/Yu6ndQEAADyUcSL47TVMahMUy/XfFkhwEIldZv2xygSAAD5KeLFcdrqGFSmqYhLPQiwUEndpv0xigQAqLoir/cp4sVx2uoYVKapiEs9CrvRMDCobvYB6WcUiZRCAEBZFLHCWrM6boy8btXiGb8TqfpBZdqKtvcbARZqLa5Ee9woUtE7KgAAmrXL1ChKv1W0i+O01TGorBsCLNRar6NIZeioAABoYL1PMdUtqKwbAizUWvMo0sTklOaazViD1Xryo6MCAJRJr5kaQFGUeUkGRS5Qe6uXje5dZLvHXVJ8NUEWpgIAyqTXIhJFLoiB+ih7lWcCLEDdVxOsY7UjAEB5RVVYO/nFo9qwZdesIKrsF7WojrJXeSZFEJkq6nRvt6l/LEwFAJRN83qfdsWaWGeMoij7kgwCLGSmyBX4eslRZ2EqAKCs2gVRSV7UFnVAFeVQ9rWDpAgiM0We7iX1DwBQB+2CqKTWGUelGq775g1a9qGrWNuFrpT9uowZLGQmi+nefkfMSP0DANRBu5mBTluXdNvHRg2oTj/levDxaUnFymBBMZX9uowAC5lJe7p30BREUv8AAFXXLohqd1HbSx/bzcApa7vQSZmvywiwkJleN/XtFYtzAQBor9PMQNxFbS99bNyAaquyFCwAekWAhcykPd3bSwoii28BAGWSZL/Vz8xAL31s1IBqlLIULAB6RYCFTKU53dttCmKRqxkCANCqCP1Wr9V2pX0DqocOD+mxJ57U9B7fe0yZChYAvaKKICqj24ozRa5mCABAqyL0W71WdVu9bFTXja3UreOv086zT9CGN7xwxmbH69csYVATlcUMFgqr13SIblMQy755HQCgXorQbw2a5l/mggVArwiwUEj9pkN0cwIv++Z1AIB6ieu35php0djmzNYSEyQB3SFFEIWUZjpE2TevAwDUS1S/JUl73Pdu5Pv+y25k816gIJjBQiF1mw7RT1Wlsm9eBwCol9Z+a46Z9rjPOIZtSYDiIMBCIXWTxjdIVSXSHAAARdNu0LC531o0tjny8awlBoqBFEEUUjdpfEWoqgSg3Mzsi2Z2r5nd1HTb4WZ2tZn9Ivz3sDzbiHpoDBpOTE7tTftb980btOxDV2nR2GatGN+6NwUwbs0wa4mBYiDAQiGtXjaq9WuWtC3pWoSqSgBK70JJr2m5bUzS99z9uZK+F/4MpGLTjgmtGN+qMy7eOWvQcPop14OPT89aZ8Va4miNz7I1IAWyRoogMtftuqlOaXxUAwQwKHe/1swWttx8kqTjwu8vknSNpPdl1ijURmuqeyeNLI3rxlZKYi1xsyJsxgw0EGAhU0meANetWjyrY2IED0ACjnL33ZLk7rvN7Mi4A81sraS1krRgwYKMmoeqiEp176SRpcFa4pnaLRvgc0LWSBFEppJcN9VNGiEApMndN7r7cndfPm/evLybg5LpJ6WdLI1oLBtAkTCDhUz1cgLsJpWQETwAKbjHzOaHs1fzJd2bd4NQTXGp7pI0Mjykx554UtN79pVjJ0sjHssGUCTMYCFT3VY+iqqmxCaKADJyhaTTwu9Pk/TtHNuCCosrVvHJNy/VzrNP0IY3vJAsjS5R+ANFwgwWMtXtuilyqQFkwcy+rqCgxRFmdpeksyWNS7rEzN4h6Q5Jb8yvhaiyThvfk6XRvU6fJZAlAixkKuoE+MrnzdOGLbt05sU7954QyaUGkAV3PyXmruMzbQhqiyAqOXyWKAoCLGSu+QQYV1Vw5KAhPfj49KzHkksNAACAIiPAQq7iUgEP2G+OhofmUoIdAABUQrf7gKL8KHKBXMVVT3poapoS7AAAoBIo3lUvzGAhN5t2TMgkecR9R48Mk0sNAAAqgeJd9cIMFnKzYcuuyODKJFIBAQBAZVC8q16YwUImovKO404qLjGaAwAAKoONkOuFGSykLi7veOSgocjjRznZAAAKbNOOCa0Y36pFY5u1Ynwr62jQERsh1wszWEgdlQIBAFURt72IRPYF4rERcr0QYCFVm3ZMtK0UeP6bl3KyAQCUBsUK0C+Kd9UHARZS0xjli0OlQABAUcXtWUSxAgCdEGAhNVGjfA2kAgIA8hYXRLVLA6RYAYBOCLCQqObOKqoEewObBgMA8tQuiIpLAzzrkhu0x33WHo69DBrGBXUAqoMAC4lp7azijIapgQAA5KXdWqq4dL89HoRVLu0NskYjgqR+ZsboF4HqIMBCYtqlBDaQGggAKIJ2a6ni0gCbNYKr68ZWzri9n5kxCmQA1cI+WEhMuwW+pqAjIjUQAFAEcWumGjNOrXsWRYnq9/qZGaNABlAtzGAhsXzwuBG/qBE+AADytG7V4llp7UNzTI8/8aTOvHinDh0e0oFDczT5+LTmmO1ND2zWCNK6WX/cbmaMAhlAtTCDVXONVIaJsFNopDL0sys9u5QDAPK2aceEVoxv1aKxzVoxvjW2P1u9bFTr1yzR6MiwTNLI8JBk0oOPT8slTU5N63fTT+n8Ny/VuW96YWz/1tqPxombGaOfBKqHGayaSzIfnF3KAQBZa549OnR4SI898aSm9wShTqciEs17Ma4Y36rJqekZ9zf6w0YWRlT/tmJ8a8f1xxa2ZcOWXTr5xaP6/s/vo58EKowAq+bi8r4nJqe0Ynxrzx0AGwcDALLSWlCiNUCSZg4atkuJ77Q+Kq5/67R+qrmk+8TklC7dPsF6ZKDiSBGsubi878Zo26BpgwAApKWb6rVSEAR1SolvV/Sinbj7R0eGNToyPCttsBHwAaguAqyai8oHb91AUaJDAAAUT7fV944eGW6bEi/1v4643eOoGgjUUy4BlpndZmY3mtlOM9uWRxvqpN2C39ZFvlGjbQ10CACAIumm+l63wU5Uf9hNKl+7x/U7Kwag3PJcg/VKd78/x9evhW52jW/NK18xvpUysgCAwohbOxVXav1pB+6nycenZxy7Ycuujn1bv+uI4x4X1T6qBgLVR5GLiuunSiAdAgCgKLoZKOymem0efRvVdYF6yivAcklXmZlL+py7b2w9wMzWSlorSQsWLMi4edXRT/53Lx1CUpsUAwAQpdNAYbezTnkFO1TXBeonrwBrhbvfbWZHSrrazH7u7tc2HxAGXRslafny5e327kMb/e4a302H0M2oIgAAgxi0UAQDgdngcwb2yaXIhbvfHf57r6TLJR2bRzvqIM1d4ztVZAIAYFBxA4JzzCKLNzXrVJodyeBzBmbKPMAys4PN7JDG95JOkHRT1u2oi36rIsVprkgYNTMmUW0QAJCcqIFCSdrj3vFinoHAbPA5AzPlkSJ4lKTLzazx+l9z9+/m0I7aSCr/uzUlMA7VBgEAg2pOOTt0eEgHDs3R5OPTmmOmPT5z5UBc8Sb2ocoGnzMwU+YzWO7+a3d/Yfj1fHf/aNZtQH+iRqhaUW0QADCo1pSzyalp/W76KZ3/5qV6yqOXZUddzLMPVTb4nIGZKNNeEVksLm03EmUSi1oBAH1r7sfazVL1UryJbUeywecMzESAVQFZVfOL69RGR4Z13djKxF4HAFAvrf1Ya3DVcPfklM5/89KuL+bZhyobfM7ATARYFdDPZsL9YIQKAJCGblLQpWCgr9eLefahygafM7APAVYF9LK4dJBUQkaoANSJmd0m6RFJeyQ96e7L821RdXVTDKF5QI+LeQBFRoBVAd3moyeRSkinBqBmXunu9+fdiKqL68fmmukpdwb0AJQKAVYFdJu6l1UqIQAAnbSWYR+aa5res2/t1fDQ3IH2bQSAvBBglVRrqt/JLx7V939+X9vUPfapAICeuKSrzMwlfc7dN+bdoCJqDZTMpMnHp9vOOrVmVExOTWtojumwg4Y6PhYAio4Aq4SiUv0u3T7RcaSvl9K2AACtcPe7zexISVeb2c/d/drmA8xsraS1krRgwYI82pirqECpoV0aelRGxfRTroP23087PnhCyq0GgHRlvtEwBtcu1a+ddasWa3ho7ozbqAIIANHc/e7w33slXS7p2IhjNrr7cndfPm/evKybmLtO1f9a+6ZNOya0Ynxr5GCfREYFgGogwCqhflP9Vi8b1fo1SzQ6MixTsH8V+e0AMJuZHWxmhzS+l3SCpJvybVXxdBMQNY5pzHbFBVcSGRUAqoEUwRJp5LlHb7/YXcdEFUAA6MpRki43MynoK7/m7t/Nt0nFE5d63nqM1Hm2i4wKAFVBgFUSrXnureiYACA57v5rSS/Mux1Z6XePxKgqts2a+6Z2s12jFLUAUCEEWCXRbuSPjgkA0K9B9khs3YC+tYrgK583Txu27NKZF+/UHDPt8dk5GKMjw7pubGXC7woA8kOAVRJxI38mddUx9Ts6CQCotkH3SIxLPW8N3KKCK7IvAFQRAVZJDFJifZDRSQBAtaW1R2Jc5sVcMz3lzmAfgMoiwCqJqDz3bkf+Bh2dBABUV9wA3hwzLRrbPCsQ6jYjIi5Ae8pdt46/Ltk3AQAFQoCVs247qtY8915G/tIanQQAlF9coYpGSl9z1oOkrjMi2NweQF0RYOWo19S9fkus08kBAOK0DuBFFaOYmt6jsy65IXIdVVxGxCCZFwBQZmw0nKN2qXtJWrdqsYaH5s64jU4OANCwetmorhtbqVvHX6enIoIoKbpIRUNURgSb2wOoK2awcpRE6l43KYaDpBcCAOqlm82Dox4Thc3tgeqjUvVsBFg5GjR1r5cUQzo5AKi3bi+COm0e3IqMCKC+qFQdjRTBHA2aupdViiEAoNwaF0ETk1Ny7bsI2rRjYtaxral9c81in5e0P6DeuBaNxgxWCrKoDChRHRAA0J1et+toznpoHaGWgsFAAisAXItGI8BKWFaVASWqAwJA3Q26J1U3F0Gs4wUQh2vRaARYCctyU19K4AJAffUyoDfoRRDreAFE4Vo0GmuwEpbGVOmmHRNaMb5Vi8Y2a8X41r0585TABYD6ihvQO+uSG2b1F1FrfofmmB5/4slZxwJAt7gWjcYMVsKSnirtNELJqCIA1FPcwF1jv6qoGa1Gmt+hw0N67Ikn9eDj07HHAkA3uBadjRmshCW9qS/VWQAAUboZuJua3qMzLt6pFeNbJWnvZsIHH7Cfpvf4rGPpWwBgcARYCUt6qrRTymFc+iAAoNqiBvTitJZlp/IXAKSHFMEUJDlV2i7lkM3dAKC+WtP+5pjtTQ+M0lxwicpfAJAeZrAKrl3KIemDAFBvq5eN7k37O/dNL+w4o9WYoUo6nR0AsA8zWAXXbv+RMy/eGfkYUjwAoDr62bw+anZK2jdDxd5WAJAeAqycddNxxqUckuIBANXW7+b1rY+TZs9QUfkLANJBimCOGh3gxOSUXLMXIXdCigcAVFu/qeDsTQMA+WEGK0ftOs5uOkFSPACg2gap9scMFQDkgwArR0mUye21A+02lx8AkD9SwQGgfEgRzFFcB5lWxzloSiIAIB1xexqSCg4A5UOAlaOsO07KugNA8bQb/GItFQCUT+VTBIucEpf1GqokUhIBAMnqtB6XtVQAUC6VDrB6LW+bhyw7TnL5AaB44ga5JiantGJ8a6EGBgEAnVU6RZCUuJnI5QeA4mk3yMVaWQAon0oHWFVNiYtbDN0JufwAUDxRg1/NWgcG++0DAADZqHSKYBVT4gZNeySXHwCKpXk9blSfJe0bGCxD6jsA1F2lZ7CqmBJH2iMAVM/qZaO6bmylRjts30EfAADFV+kAq4opcVVNewQAdB4YpA8AgOKrdIqgVL2UuCqmPQIAAp2276APAIDiq3yAVTXrVi2ekX8vlT/tEQCKyMxeI+lTkuZK+oK7j2fxuu0GBukDAKD4CLBKJuvNiQGgjsxsrqTPSHq1pLsk/djMrnD3n+XZLvoAACg+AqwSqlraIwAU0LGSfunuv5YkM/uGpJMk5RpgSfQBAFB0lS5yAQBAn0Yl3dn0813hbTOY2Voz22Zm2+67777MGgcAKC4CLAAAZrOI23zWDe4b3X25uy+fN29eBs0CABQdARYAALPdJemYpp+fKenunNoCACgRAiwAAGb7saTnmtkiM9tf0lskXZFzmwAAJZBLgGVmrzGzXWb2SzMby6MNAADEcfcnJb1b0hZJt0i6xN1vzrdVAIAyyLyKYFFL3wIA0MzdvyPpO3m3AwBQLnnMYO0tfevuT0hqlL4FAAAAgFLLI8Ci9C0AAACASsojwKL0LQAAAIBKyiPAovQtAAAAgEoy91mTR+m+oNl+kv63pOMlTSgohfsX7aozmdl9km7PpoUzHCHp/hxetyjq/P7r/N6ler//Or93Kd/3/yx3L23KQlNfVfe/oW7wGbXH59MZn1F7fD6d9fMZddVPZV5F0N2fNLNG6du5kr7YqfRtXh2umW1z9+V5vHYR1Pn91/m9S/V+/3V+7xLvfxCNvorPsDM+o/b4fDrjM2qPz6ezND+jzAMsidK3AAAAAKopl42GAQAAAKCKCLDa25h3A3JW5/df5/cu1fv91/m9S7z/JPAZdsZn1B6fT2d8Ru3x+XSW2meUeZELAAAAAKgqZrAAAAAAICEEWAAAAACQEAKsGGb2GjPbZWa/NLOxvNuTJTO7zcxuNLOdZrYt7/akzcy+aGb3mtlNTbcdbmZXm9kvwn8Py7ONaYp5/+eY2UT4N7DTzE7Ms41pMbNjzNiyuqIAAAT3SURBVOz7ZnaLmd1sZu8Jb6/877/Ne6/F7z4tde47ulG3/qUbde+DOqlzH9WtOvdl3cijv2MNVgQzm6tgM+RXS7pLwWbIp7j7z3JtWEbM7DZJy929FhvUmdnLJT0q6Uvu/oLwto9LesDdx8OLpMPc/X15tjMtMe//HEmPuvsn8mxb2sxsvqT57v4TMztE0nZJqyWdror//tu89zepBr/7NNS97+hG3fqXbtS9D+qkzn1Ut+rcl3Ujj/6OGaxox0r6pbv/2t2fkPQNSSfl3CakxN2vlfRAy80nSboo/P4iBf8RKynm/deCu+9295+E3z8i6RZJo6rB77/Ne0f/6DvQs7r3QZ3UuY/qVp37sm7k0d8RYEUblXRn0893qV4XHi7pKjPbbmZr825MTo5y991S8B9T0pE5tycP7zazn4bpGZVPKzCzhZKWSbpeNfv9t7x3qWa/+wTVve/oBv1Ld2p1DuoT56kIde7LupFVf0eAFc0ibqtTLuUKd3+RpNdK+rtweh718llJz5a0VNJuSefm25x0mdnTJF0q6Qx3fzjv9mQp4r3X6nefsLr3Hd2gf0ESOE9FqHNf1o0s+zsCrGh3STqm6ednSro7p7Zkzt3vDv+9V9LlCtJe6uaeMGe3kbt7b87tyZS73+Pue9z9KUmfV4X/BsxsSMEJ96vufll4cy1+/1HvvU6/+xTUuu/oBv1L12pxDuoX56nZ6tyXdSPr/o4AK9qPJT3XzBaZ2f6S3iLpipzblAkzOzhcACgzO1jSCZJuav+oSrpC0mnh96dJ+naObclc44Qcer0q+jdgZibpAkm3uPt5TXdV/vcf997r8rtPSW37jm7Qv/Sk8uegQXCemqnOfVk38ujvqCIYIyzV+ElJcyV90d0/mnOTMmFmf6hgVFGS9pP0taq/dzP7uqTjJB0h6R5JZ0vaJOkSSQsk3SHpje5eyUW2Me//OAVT5i7pNknvauRxV4mZ/amk/ynpRklPhTd/QEFudqV//23e+ymqwe8+LXXtO7pRx/6lG3Xvgzqpcx/VrTr3Zd3Io78jwAIAAACAhJAiCAAAAAAJIcACAAAAgIQQYAEAAABAQgiwAAAAACAhBFgAAAAAkBACLCBHZnaMmd1qZoeHPx8W/vwsM/uumU2a2ZV5txMAUE9t+qlXmNkPzexmM/upmb0577YCRUGZdiBnZvaPkp7j7mvN7HOSbnP39WZ2vKSDFOzL8Gf5thIAUFdR/ZSkSyW5u//CzI6WtF3Sf3L3yRybChQCARaQMzMbUtAxfVHSOyUtc/cnwvuOk/ReAiwAQF7a9VNNx9wg6Q3u/oscmggUyn55NwCoO3efNrN1kr4r6YTWTgsAgDx16qfM7FhJ+0v6VR7tA4qGNVhAMbxW0m5JL8i7IQAARIjsp8xsvqQvS3q7uz+VR8OAoiHAAnJmZkslvVrSSyWdGXZWAAAUQlw/ZWZPl7RZ0j+5+3/k2ESgUAiwgByZmUn6rKQz3P0OSRskfSLfVgEAEIjrp8xsf0mXS/qSu38zzzYCRUOABeTrnZLucPerw5//u6TnheVv/6ekb0o63szuMrNVubUSAFBXkf2UpPdLermk081sZ/i1NK9GAkVCFUEAAAAASAgzWAAAAACQEAIsAAAAAEgIARYAAAAAJIQACwAAAAASQoAFAAAAAAkhwAIAAACAhBBgAQAAAEBC/n+4L/+gRokM4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x396 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1,ax2) = plt.subplots(1,2 , figsize = (12,5.5) )\n",
    "\n",
    "ax1.plot(x1, y, 'o')\n",
    "ax1.set_xlabel('X1')\n",
    "ax1.set_ylabel('Y')\n",
    "ax1.set_title('Contant variance of residuals')\n",
    "\n",
    "ax2.plot(x2, y2, 'o')\n",
    "ax2.set_xlabel('X2')\n",
    "ax2.set_ylabel('Y')\n",
    "ax2.set_title('Non contant variance of residuals')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr.fit(x1.reshape(-1,1),y.reshape(-1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr2 = LinearRegression()\n",
    "lr2.fit(x2.reshape(-1,1),y.reshape(-1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGECAYAAAAr/IHVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XecVNX5x/HPs8sCi4W1oMIqQtRAVAwowRiMwYLYXTHGXtKwxhJF0FjWgqDY8lNjr7Fh1KxGQSyIhVgCgqAiISqoCyoqS3N32XJ+f9y7O2VnZmd2p8/3/Xrti7ll7j07wD3znPOcc8w5h4iIiIiIiHReUaYLICIiIiIiki8UYImIiIiIiCSJAiwREREREZEkUYAlIiIiIiKSJAqwREREREREkkQBloiIiIiISJIowBJJMjO72MzuyXQ54mVmw81ssZmtNbOKFN6nr3+P4ijHK83s4STdy5nZ9sm4loiIeMxsmpmdnOlyxMvMTjezr/26Z7MU3ud4M3sxxvGZZvaHJNxnhJl92dnrSOopwJIOM7PjzGy2/+Ba7j9490zCdR8ws6uTUUb/ekl5sMXLOXeNcy5t90uCK4FbnXMbOueqUnUT59zn/j2aUnUPEclfZrbE/7K8QdC+P5jZzAwWKyH+77BfEq+X1sYk59yBzrkH03W/zjCzEuBGYH+/7vkuVfdyzj3inNs/VdeX3KMASzrEzP4M3AxcA2wJ9AX+BhyeyXJlmpl1yXQZOmBb4MN4TszR309E8kcX4JxMF6LQmCfXvjNuCXRH9ZtkQK79Z5EsYGY98Xo9znTOPe2cW+eca3DO/cs5N9Y/p5uZ3Wxmy/yfm82sm39shJl9aWbnm9k3fu/Xb/1jY4DjgQv9nrF/+fvHm9knZrbGzD4ysyOCynOKmb1pZteb2Uoz+8zMDvSPTQB+CdzqX+/WCL/PC2Z2Vti+981stP/6r2b2hZmtNrM5ZvbLoPMqzexJM3vYzFYDp4SnupnZP8zsKzNbZWavm9lOQcceMLPbzOx5/3d7x8y2Czq+k5m9ZGbf+y23F/v7i4I+k+/M7Akz2zTG39kfzex//nWeNbM+/v5PgB8B//I/n24R3rvEzMaZ2XxgnZl1MbM+ZvaUma3wP++zg84f5vdsrvbLfKO/v5/f2trF3+5vZq/5v/dLwOZB12iTBhHc8uvf4y0zq/H//dxqZl2j/O4H+f9m1phZtZldEO1zEpGsNxm4wMzKIh00s1+Y2X/85+1/zOwXQcdmmtlVZjbLfx68aGabR7qOf/7hZjbPf5Z9YmYH+Pv7+M/R7/3n6h+D3lPpP48f8u/xoZkN9Y/9Ha8xsuV5e6G/v0N1hJm97p/2vn+9o8PK381/Ru4ctK+XmdWa2RZmtomZPec/x1f6r7cO+7wmmNks4AfgRxaUEWJm25nZDL8O+tbMHgn+e/Gf2ReY2Xz/d5tiZt3j+Hx7mtm9/rO92syutuip5RG/a5jZj4FF/mk1ZjYjwntb6qTfm9nnwAx//8/N7N/+Z/e+mY0Ies8pZvap/3fxmZkdH7T/zaDzRprZx/7vfStgYf9GHo5Qjpa68bdmttC/x6dmdmqk390/d5z/Ga0xs0Vmtm+0cyXNnHP60U9CP8ABQCPQJcY5VwJvA1sAvYB/A1f5x0b4778SKAEOwnt4b+IffwC4Oux6RwF98BoFjgbWAb39Y6cADcAfgWLgdGAZYP7xmcAfYpT1JGBW0PaOQA3Qzd8+AdgMr+X0fOAroLt/rNK/d4VftlJ/38NB1/sdsBHQDa/Xb17QsQeA74Fh/vUfAR73j20ELPfv2d3f3t0/dq7/+W7tX/dO4LEov98+wLfArv65twCvBx1fAuwX4/NZAswDtvF/vyJgDnAZ0BUvQPsUGOWf/xZwov96Q+Dn/ut+gGv5d+Ofd6Nfpr2ANS2fm/9v5MsI5djPf70b8HP/M+sHLATODTrXAdv7r5cDv/RfbwLsmun/Q/rRj34S/2l5BgBP49cRwB+Amf7rTYGVwIn+s+FYf3sz//hM4BPgx/6zbCYwKcq9hgGrgJH+M68cGOgfew0vY6M7MBhYAezrH6sE6vDqtWJgIvB2+O8Qdq8O1RH+8dZnXZTf4z5gQtD2mcAL/uvNgCOBHv79/wFUBZ07E/gc2Mm/dwlB9Smwvf/5dMOr518Hbg77Xd/Fq7s39Z/Tp8Xx+Vbh1Wkb4H2HeBc4NcrvF+u7Rj+C6pwI7205/pB/r1K/HN/5f39Ffvm+86+9AbAaGOC/vzewk//6FOBN//Xm/nm/9j+z8/C+87R8bpWEfkcIKSdwMLAdXlD2K7zvR7v6x0bg143AAOALoE/QdbbL9P9T/fh/r5kugH5y7wevh+mrds75BDgoaHsUsMR/PQKoDX7oAd8Q+CL+AGEBVoTrzwMO91+fAvwv6FgP/2G1lb89k9gB1kZ4Adu2/vYE4L4Y568Efuq/riQoWAna93CU95b5ZesZ9LveE3T8IOBj//WxwNwo11mIX6H7273xAr02FQlwL3Bd0PaG/rn9/O0ltB9g/S5oe3fg87BzLgLu91+/DlwBbB52TmslgteK2whsEHT8UeIMsCKU8Vzgn0HbwQHW58CpwMaZ/r+jH/3op+M/BAKsnfG+nPciNMA6EXg37D1vAaf4r2cClwQdOwM/2IhwrzuBmyLs3wZoAjYK2jcReMB/XQm8HHRsR6A2/HeI8TvGXUf42+0FWPsBnwZtzwJOinLuYGBl0PZM4Mqwc2YSpT7Fa2icG7S9BDghaPs64I52Pt8tgXqgNGjfscCrUe4Z67tGP+ILsH4UtG8c8Pew86YDJ+MFWDV4QWlp2DmnEAiwTiI0qDbgS+IMsCKUswo4x389gkCAtT3ed6f9gJJ0/B/UT/w/ShGUjvgO2Nxi5yv3AZYGbS/197VewznXGLT9A94X/4jM7CQ/laDGzGrwKtjg1I6vWl44537wX0a9XjDn3BrgeeAYf9cxeK2ELfc+3++uX+Xfu2fYvb+IUe5iM5vkpz+sxqtwiFZ2Qj+HbfAqj0i2Bf4Z9HksxKv0t4xwbsjfhXNuLd7fYXm0ckcQ/DtuC/Rpubd//4uD7v17vBbij81L0TkkSplWOufWBe1bGuG8iMzsx346y1f+53oNoZ9psCPxvpQsNS8lcY947yMi2cc59wHwHDA+7FB4vYO/Hfysi/a8DRft+dsH+N6vN+K9R/do9WUn64h4zABKzWx3M9sWL4j6p3/vHmZ2p5kt9e/9OlAWlo4Xq37bwswe91PUVgMP0/Y5nGj9ti1er8/yoPrlTrweqkja+64Rj/D67aiw+m1PvIyZdXgZNKf55XvezAZGKVPrNZ0XDUX9HMOZ2YFm9rZ5Kag1ePVXm/rNOfc/vMbFSuAb/+8i0d9dUkQBlnTEW3gpELGm9F6G96Bq0dffFw8XvOFXCncDZ+GlepQBHxCU05zI9aJ4DDjW//JdCrzq3/uXeC1av8FLYSzDazkNvnes6x+HN/HHfniBWT9/fzxl/wIvTSDasQOdc2VBP92dc9URzg35uzBvBq7NgEjnRhP8O34BfBZ2742ccwcBOOcWO+eOxasQrwWetKBZv3zLgU3C9vcNer0OryeypczFeK3VLW4HPgZ2cM5tjBfgRfxMnXP/cc4d7penCngi7t9aRLLV5Xhp4cGBTXi9A95zJZFnXYtoz99lwKZmtlEH7xFeX3Smjmj/Zs414z3zjvXv9VxQcHg+XprZ7v5zdK8I945Vv030j+/iv/+EBMod7fP9Aq8Ha/Og+mVj59xOEc6Fzn3XaBFev/09rH7bwDk3CcA5N905NxIva+RjvO8m4ZbjBZCAN0FI8DZh9RuwVdC53YCngOuBLf3vHFOJXr896pzbE+8zcHh1rmQBBViSMOfcKrzxN7eZWYXfClbit7pc55/2GHCJeQNqN/fPj3eNo6/xxvW02ADvwbECvAGgeD1Y8Qq/XiRT8R5QVwJT/EoJvPTBRv/eXczsMmDjBO69EV5l8R3eA/WaBN77HLCVmZ3rD9rdyMx294/dAUzwg8+WgcvRZnB8FPitmQ32H97XAO8455YkUJZg7wKr/cG1pX4L7M5m9jO/LCeYWS//M6zx3xMyNbtzbikwG7jCzLqaN73/oUGn/Bev1fdg86bavQQvz7/FRng57mv9FsTTIxXUv/bxZtbTOdfgv0fTxIvkOL/1fgpwdtDuqcCPzVtCpIt5kz7siPcsTdS9eM/Nfc2bVKjczAY6577AG+cz0cy6m9kueL32j8S8WkB4fdSZOiLS9SJ5FK/n5Xj/dfC9a/EmgdgUL2hNxEbAWv/95cDYBN4b7fNdDrwI3GBmG/vHtjOzX0W5Tme+a0TyMHComY3y67bu5k26tLWZbWlmh/kNg/V4v3uk+uR5YCczG+33XJ5NUBCFN8RhL/PWhuyJl2LfoiteXbcCaDRvwq6I07+b2QAz28ev1+vw/i5Vv2UJBVjSIc65G4E/433xXYHX6nMWXg8BwNV4X6DnAwuA9/x98bgX2NHvnq9yzn0E3IDXc/Y1MAgvjzxefwV+bd4sSf8X5fepxxs4vR+hFdB0YBreF/6leA+xuLv68QbPLsVr3fwIbzBuXPxWxpF4gcdXwGJgb//wX4FngRfNbI1/3d2jXOcV4FK8VrHleK2Gx0Q6N85yNfllGgx8hjeBxj14ra/gTYLyoZmt9ct5jHOuLsKljvPL/D1exf5Q0D1W4Y2PuAfvs1uHl8Pe4gL//WvwWhCnxCjyicASP4XlNLxWVhHJfVfiNcAB4Lx1jg7B65n5DrgQOMQ5922iF3bOvQv8FrgJL2vhNQI9Jcfi9TQtw0u3u9w591Kcl56IFxDUmDejaYfrCF8l8KB/vd9E+V3ewXuG9sGrz1rcjJex8a1/3xcSvPcVeJMnrcILKp6O943tfL4n4QUaH+GNeX4Sr8coks5814hUri/wehQvJvDdZize9+UivH9by/DqrV/h1VPh1/gWb2KuSXj/Dncg6DuL/29lil/mOQQ1APj1/tl4vY4r8eq5Z6MUt5t/j2/xviNs4ZdbskDLLGsiIiIiIiLSSerBEhERERERSRIFWCIiIiIiIkmiAEtERERERCRJFGCJiIiIiIgkSayFYrPG5ptv7vr165fpYoiISArNmTPnW+dcr/bPzE6qq0RE8lu89VROBFj9+vVj9uzZmS6GiIikkJktzXQZOkN1lYhIfou3nlKKoIiIiIiISJIowBIREREREUkSBVgiIiIiIiJJogBLREREREQkSRRgiYiIiIiIJIkCLBERERERkSRRgCUiIiIiIpIkCrBERERERESSRAGWiIiIiIhIkijAEhERERERSZIumS6AiIiIiIgUnqq51UyevohlNbX0KStl7KgBVAwpz3SxOk0BloiIiIiIpFXV3GouenoBtQ1NAFTX1HLR0wsAcj7IUoqgiIh03PL3YfFLmS6FiIjkmMnTF7UGVy1qG5qYPH1RhkqUPOrBEhGRjnn3bph6gff68howy2x5REQkZyyrqU1ofy5RD5aIiCSu6sxAcPWbv+dscGVm3c3sXTN738w+NLMr/P0PmNlnZjbP/xmc6bKKiOSTPmWlCe3PJerBEhGRxNy8C9Qs9V6f/m/YcqfMlqdz6oF9nHNrzawEeNPMpvnHxjrnnsxg2URE8tbYUQNCxmABlJYUM3bUgAyWKjkUYImISHwa18PVvQLbF34GPTbNXHmSwDnngLX+Zon/4zJXIhGRwtAykYVmERQRkcK0dgVcv31g+9Jvobgkc+VJIjMrBuYA2wO3OefeMbPTgQlmdhnwCjDeOVefyXKKiOSbiiHleRFQhdMYLBERiW35+4HgqtdAqFyVN8EVgHOuyTk3GNgaGGZmOwMXAQOBnwGbAuMivdfMxpjZbDObvWLFirSVWUREspcCLBERie6Dp+DOvbzXQ38HZ76T2fKkkHOuBpgJHOCcW+489cD9wLAo77nLOTfUOTe0V69ekU4REZECk7IAy8y2MbNXzWyhPzPTOf7+SjOrDpqZ6aBUlUFERDrhpcvhyd95rw+7BQ65KbPlSQEz62VmZf7rUmA/4GMz6+3vM6AC+CBzpRQRkVySyjFYjcD5zrn3zGwjYI6ZtaxGeZNz7voU3ltERDrjvgPg87e81797EfruntnypE5v4EF/HFYR8IRz7jkzm2FmvQAD5gGnZbKQIiKSO1IWYDnnlgPL/ddrzGwhkH+j2ERE8klzM1y5SWD7vI+gZ/4+up1z84EhEfbvk4HiiIhIHkjLGCwz64dXgbUk759lZvPN7D4z2yTKezRwWEQkDlVzqxk+aQb9xz/P8EkzqJpb3bEL1a0ODa7+8lVeB1ciIiKpkPIAy8w2BJ4CznXOrQZuB7YDBuP1cN0Q6X0aOCwi0r6qudVc9PQCqmtqcUB1TS0XPb0g8SBr8UswaRvvdfeecHkNlJQmvbwiIiL5LqUBlpmV4AVXjzjnngZwzn3tT4nbDNxNlJmZRESkfZOnL6K2oSlkX21DE5OnL4r/ItP/Ao/82ntd0gPGfw5mSSyliIhI4UjZGCx/5qV7gYXOuRuD9vf2x2cBHIFmZhIR6bBlNbUJ7W/jhoGwxn8kl24C45bEfe+qudVMnr6IZTW19CkrZeyoAXm5YKSIiEgiUjmL4HDgRGCBmc3z910MHGtmgwEHLAFOTWEZRETyWp+yUqojBFN9yuJI76vsGXg9bAwcNDnu+7akJrb0nrWkJgIKskRE8owa1BKTylkE38Sb3jbc1FTdU0Sk0IwdNSAk0AEoLSlm7KgB0d/kHFxRFtg+5jEYmNiShLFSE1XpiojkDzWoJS6VPVgiIpIiwa2JPUtL6F5SRM0PDe23LNbWwLXbBrbPXQBlfRO+f6dTE0VEJCeoQS1xCrBERHJMeGtiTW0DpSXF3HT04NiVXfUcuDtoeadLVkCXrh0qQ6dSE0VEJGeoQS1xaVkHS0REEhNrbasOzRz49u2hwVXlqg4HV+ClJpaWFIfsazc1UUREck60hjM1qEWnAEtEJMu0t7ZVwq2J9x8ML4wPbFeu6nQZK4aUM3H0IMrLSjGgvKyUiaMHKV1ERCTPqEEtcUoRFBHJMu3luyeUnhc8U+AO+8Px/0haOSuGlCugEhHJcy3Pec0iGD8FWCIiWaa9HqpIMweWFBvr6hvpP/75QOX3zI6BNx98I/zs9yktt4iI5Cc1qCVGAZaISJZpr4cqvDWxrEcJa+saqaltAGBFzerQ4Oq0N2GrQakvuIiIiGgMlohItokn371iSDmzxu/DZ5MOpkfXLjQ0OwC2ta/4b/eTA2+86EsFVyIiImmkAEtEJMskOoFES+rgQUVv81q3P7fu71f3CFUfrU5HkUVERMSnFEERkSyUSL57n7JSzl77V47uMrN1X7+6RwG46OkFrdcTERGR1FMPlohIjntmg6tbg6tPmnu3BlcQx/pYIiIiklTqwRIRyWWVPdncf3lLYwU3NP6mzSlR18cKUjW3WlPwioiIJIECLBGRXNTcDFduEtgeM5PHH1oJ8a6PFaRlYeOWad9bFjYGpRaKiIgkSimCIiK5prYmNLi68DPoMyTi7IPB62MNnzSDqrnVbS4Xa2FjERERSYx6sEREcslXC+COPQPbl62EIq+trL31saL1TLW3sLGIiIjETz1YIiK5Yt5jocFV5arW4KpFtPWxWkTqmYqWQtheaqGIiIi0pQBLRCQXPHMmVJ0GwLvNAxne/Z8R0/2CxdszFc/CxiIiIhIfpQiKiGS7a/tD7fcATG74Dbc1VcD69iei6FNWSnUck16EpxZqFkEREZGOU4AlIpLNKnu2vjx+/UXMah7Uut2S7hctEBo7akDI7IAQvWcqkYWNRUREJDoFWCIiWaDNOlT7b0/Fs4Fg6hd1t7CMzdq8L9ZEFOqZEhERST8FWCIiGRa+DtX3NSvZ+JlTwPwTLvkGu35Wh9a4Us+UiIhIemmSCxGRDAteh6oXK5nS9Sp+xXvc0OWP3kyBXbppIgoREZEcoR4sEZEkaZPmF2c6Xkua3znFT3FeyVPUuxL+2HA+r9bvyvn+OUr3ExERyQ0KsEREkiA8zS/aor6R9Ckr5W8/nM9Piz4F4Jj1lzDX7UB5hNn+FFCJiIhkN6UIiogkQXCaX4tIi/pGMqvuiNbg6uz1ZzHX7UBJsSn9T0REJAepB0tEJAniXdQ3RHMTXLlp6+bI+utY7Lb2NlwySyciIiLpoh4sEZEkiDabX9RZ/tZ9GxJc7Vx3TyC4AhqaXVy9XyIiIpJd1IMlItIB4RNa7D2wF0/NqY5rUV++nA337Nu62b/uEVzrnOwBMXu/REREJCupB0tEJEEtE1pU19Ti8Ca0eGpONUfuVk55WSkGlJeVMnH0oLaTUrx7d0hwReUq+pT1iHif9ta4ks4zs+5m9q6ZvW9mH5rZFf7+/mb2jpktNrMpZtY102UVEZHcoB4sEZEERZvQ4tWPVzBr/D7R3zjlBFj4L+/1jw+E4x4HYOyoASEzEILWuEqjemAf59xaMysB3jSzacCfgZucc4+b2R3A74HbM1lQERHJDerBEhFJULTUveqaWoZPmkHV3Oq2Byt7BoKrAya1BlfgTb8+cfSg9nu/JOmcZ62/WeL/OGAf4El//4NARQaKJyIiOUg9WCIiCepTVkp1jCArZP0r5+CKssAJv5sOfX/e5n1a4ypzzKwYmANsD9wGfALUOOca/VO+BCL+5ZjZGGAMQN++fVNfWBERyXrqwRIRSdDYUQMoLSmOerx1/at134UGVxcsjhhcSWY555qcc4OBrYFhwE8inRblvXc554Y654b26tUrlcUUEZEcoR4sEZEowmcKHDtqQEhP0+Tpi6L2ZG27ejZMPiKw49LvoFiP3GzmnKsxs5nAz4EyM+vi92JtDSzLaOFERCRnqAdLRCSCSDMFXvT0gtbxVRVDypk1fh/KI8z0N77LozzadUJgR+UqBVdZysx6mVmZ/7oU2A9YCLwK/No/7WTgmcyUUEREco0CLBGRCKLNFBi++G94uuCS7sdxWpfnWreHd/8n/cc/H33yC8m03sCrZjYf+A/wknPuOWAc8Gcz+x+wGXBvBssoIiI5RE2qIiIRRJspMHx/cLrgrLpASuAPpb3Zbe3N1Prnt5n8QrKCc24+MCTC/k/xxmOJiIgkRD1YIiIRRFvkN9L+iiHlIcEVw89lpPtbXD1gIiIikl8UYIlIQauaW83wSTPapPFFmikw4uK/dau8Na5a/G46jLwi7h4wERERyS9KERSRgtUykUVLT1OkNL5Iswi2+ngqPH5sYPu8D6Hn1kD0tbKi9YyJiIhIflCAJSIFK9ZEFi3TsUcdLzX1Qnj3zsD25TVg1ro5dtSAkOANovSAiYiISF5RgCUiBavDaXzBKYHgTcNO23WzjtytnFc/XhG9B0xERCRB0dZoLCTZ/hkowBKRgtWhNL4YwVV4uuFTc6qZOHpQVj30RUQkd8WT2p7vcuEz0CQXIlKw4p7IokVwcLXhlq3BFcS/bpaIiEhHqa7Jjc9APVgiUrBiTWQRnH6wdc+uvFH/68AbR14Jw88JuZZmDRQRkVRTXZMbn4ECLBEpaJEmsghOP+hny5lZf37g4BnvwBYD21xHswaKiEiqqa7Jjc8gZSmCZraNmb1qZgvN7EMzO8ffv6mZvWRmi/0/N0lVGUSksERb0ypRLekH+xbNYWa3QHC1V7d/RAyuoAPphiIiIglSXZMbn0Eqe7AagfOdc++Z2UbAHDN7CTgFeMU5N8nMxgPjgXEpLIeIFIBkDnpdVlPLxC53c2yXV1v39at7FKtriPqeuNbNEhER6QTVNbnxGaQswHLOLQeW+6/XmNlCoBw4HBjhn/YgMBMFWCLSSe2taZWIz7ofF7Ldr+5RoP30g5jrZomIiCSB6prs/wzSMgbLzPoBQ4B3gC394Avn3HIz2yLKe8YAYwD69u2bjmKKSA5L2qDXoJkC5zf357D1E4DY6QfZvh6HiIhIuqhOTMM07Wa2IfAUcK5zbnW873PO3eWcG+qcG9qrV6/UFVBE8kK03qWEBr0GBVcLdh7H6T1uxIDystKo61m1pCZW19TiCKQmdnT8l4iISK5SnehJaQ+WmZXgBVePOOee9nd/bWa9/d6r3sA3qSyDiBSGsaMGhIzBggQGvTbUwoStAtunvsGg3rswK477JjM1UUREJJepTvSkLMAyMwPuBRY6524MOvQscDIwyf/zmVSVQUQKRyKDXoPTF3bbeDVPrj8tcPDiZdB1g7jvmwvrcYiIiKSD6kRPKnuwhgMnAgvMbJ6/72K8wOoJM/s98DlwVArLICIFJJ5Br8GzDY4omscD668LHKxclfA9c2E9DhERkXRQnehJ2Rgs59ybzjlzzu3inBvs/0x1zn3nnNvXObeD/+f3qSqDiEi4lvSF87r8gwe6esHVcrcpw7v/s0PXy4X1OERERNJBdaInLbMIioikU6wZjJbV1PJM10v4adGnAPy9cT8ubfwdVt+x9IVcWI9DREQkHVQnehRgiUheaW/B4eA1rs5efxbPNv8C6Fz6QravxyEiIpIuqhMVYIlIHqmaW835T7xPk3Mh+2sbmrjhhY+oeGbH1n0j669jsdsaKMz0BREREUmNlK+DJSKSDi09V+HBFcCmrOaN+l+3bj930Lv80HOHdte4EhEREUmUerBEJC9EWnsDYLD9j6pulwV2XF7DIWYcMkw9ViIiIpJ8CrBEJGvFmqwiXKQ1Nk4ofomrS+4P7OjANOwiIiIiiVCAJSJZqb3JKsKFr71xW8nNHFz8LgBfbTmCrU7XmuYiIiLZKpFG1WynMVgikpUipfzVNjQxefqiiOcHr73xSbfjW4Or+TtfpOBKREQki7U0qlbX1OIINKpWza3OdNE6RD1YIpKVIqX8ARFXiAe/V8s5Kp7dqXXf63s+zF77HZqS8omIiEhyxGpUzcVeLPVgiUhWirYulUHkFq2mRiqqrw9sn/9fBVciIiI5IFqjarT92U4BlohkpbGjBmAR9jtomyZYvwYeOxrm3A97ngeXrYSNtkxHMUVERKSTojWqRtuf7RRgiUhWqhhSTtsVrTzLamqpmlvN8Ekz2GP8QyyetCfNn7wKh/4V9quEIj3aRERL3sKeAAAgAElEQVREckXwOOoWpSXFjB2Vm0uqaAyWiGSt8rCZAVuU9SjhoqcXcFjzy1zb/W7WNJcypvlCDikaSUUGyikiIpJP0j2jX8u182UWQQVYIpIUqXgYjx01IGSqdvBatJyDhcVHg9/YddT6y/nY9WVhjg6GFRERyRaJLpOSLBVDyvOmDlcejYh0WqqmV60YUs7E0YMoLyvF8Hq0Jo4exDx3VOs5NzeO5mPXF8jdwbAiIiLZItFlUqQt9WCJSKelcnrVkBat+rUwMXC9X9dfxmw3sHU7VwfDioiIZIt8m9EvE9SDJSKdlpaH8dJ/hwRXuzXcHRJclRRbzg6Glcwxs23M7FUzW2hmH5rZOf7+SjOrNrN5/s9BmS6riEg65NuMfpmgAEtEOi3lD+MHD4P7D2zdrDrsQ1a7DUPPiTbloEhsjcD5zrmfAD8HzjSzHf1jNznnBvs/UzNXRBGR9Mm3Gf0yQQGWiHRaSh/GlT3hs9eCtlcx+cX/0tAcGlE1NDvlh0vCnHPLnXPv+a/XAAuB/BhlLSLSAdHGP2dyAoqWpVn6j3+e4ZNmdHqMd6ppDJaIdFo806t2aJbByp5h26sA5YdLaphZP2AI8A4wHDjLzE4CZuP1cq2M8J4xwBiAvn37pq2sIiKplE0z+mVqVsPOUIAlIkkR62HcoYdjcHC1wRYwdnHrZp8o62MpP1w6ysw2BJ4CznXOrTaz24Gr8JJPrwJuAH4X/j7n3F3AXQBDhw5VoqqISJJ1eiKt2pUwcxJstj0M+2OKShlKKYIiknIJTfna3BQaXO13RUhwBcoPl+QysxK84OoR59zTAM65r51zTc65ZuBuYFgmyygiUqg6nLXS3AxzHoRbdoN374LV6UsrVA+WiCRdeDpgpN4miPBw/O4TuGXXwPYZ78AWAwmXbyu+S+aYmQH3AgudczcG7e/tnFvubx4BfJCJ8omI5LIODQ8I06GslS9nw9QLYNlc6LsHHHgt9P5posXvMAVYIpJUkdIBjciT/IU8HN+5C6aNDWxf+i0Ul0S9Tzblh0tOGw6cCCwws3n+vouBY81sMN4/3SXAqZkpnogUumQEKZmQrLFTY0cNCLkOxMhaWfM1vFwJ7z8KG24Fo++GQUeBWad+l0QpwBKRpIqUDuigTZAV8nC8ZSh8F5QG6E9mIZJqzrk38f55htO07CKScbk4wUOLTo+d8sWVtdLUAO/c6Y21aqyD4efCXhdAt42S8rskSgGWiCRVtJzo4OBqkx4lXH7oTt7DMcpMgSIiIoUuWUFKJiRzxt+YWSufzIBp4+Db/8L2I+GASbD59gnfI5kUYIlIUsUac9WirqHZe6HgSkREJKpcXpYk5TP+rlwK0y+Gj5+DTfrDsVNgwAHJuXYnaRZBEUmqSDP8hattaKLimR0DOwYfHzO4yrUFBkVERJIhWjCSC8uSpGzG34ZaeHUi3DbM673a51I44+2sCa5APVgikmThudLhk1uUsYZ53YPmCzixCrbbO+r1cjn/XEREpDMSmuAhyyR9xl/nYOG/YPpfYNXnsNNo2P8q6Ll1EkudHAqwRCTpgnOlh0+a0ZoisHfRXO7vOjlw4rilUFoW81q5nH8uIiLSGbm+LEnSZvz95mN4YRx8OhO22AlOfg76/7Lz100RBVgiklItrW/XcROHFr/duv+SwW9ydTvBFeR2/rmIiEhnFfSyJHWrYOa18O6d0HUDOPA6GPp7KM7uECa7SyciOa9iSHnoeCugX92jlM6pZui2m7ZbaaR8kKyIiIhkl+Zmby2rlyth3bew28neWKsNNs90yeKiSS5EJLXCZgrsV/coEEjza0/KBsmKiIhI9qmeA/eOhGfOhE36wZhX4dC/5kxwBerBEpFOirnCfFBwtaC5H4euvybkvfGk+eV6/rmIiIjEYe0KeOUKmPswbNALKu6AXY6GotzrD1KAJSIdFm2GP2tu4PB//bT1vIldzuTOtcPbvD/eNL+Czj8XERHJZ00N8J97vKnXG9bBHmfCr8ZB940zXbIOU4AlIh0WaYa/zRuXcfi/jg7sOGc+P1nShdIcnWZWREREUuSz12HqhbBiIWy3DxxwLfT6caZL1WkKsESkw8JT/A4vepO/dv1bYMdlK6GoiIpNvM3OpvnFTEcUERGR3FDzBbx4CXxUBWV94ehHYODBYJbpkiWFAiwR6bDgGf6OL36ZCSX3BQ5Wrgo5t7NpflpwWEREJMc11MG/b4E3bgAcjLgYhp8NJfk1M3DujRoTkawxdtQASoqNM4urQoKrqsM/Svq9Yi04LCIiIlnMOfj4ebhtGLx6NewwEs76D4wYl3fBFagHS0Q66eMux1FsDoDKhpN4xB3I5BTcRwsOi4hILlFau+/bxTBtHHzyCvQaCCc9Az8akelSpZQCLBHpGOe8BYT9dOl96q/nU9cHcEyevijplYgWHBYRkVyhtHagbjW8fh28fTuU9IBRE2HYH6G4JNMlSzmlCIpIG1Vzqxk+aQb9xz/P8EkzqJpbHXpC/Vq4oqx1c2jd7X5w5UlFr5IWHBYRkVxR0GntzsH7j8OtQ73xVj89Bv70HuxxRkEEV6AeLBEJ026r27eLvYemb7u6v9NEaOCTil4lLTgsIiLZJloaYMGmtS+bB9MuhC/egT67wjGPwda7ZbpUaacAS0RCxGp1q+g2B544sXV/1eEf0TWN61tpwWEREckWsRokCy6tfd13MOMqmPMA9NgMDrsVBh8PRYWZLKcAS0RCRKoQAH679i54Ypq3seXOcPosKvxj6lUSEZFc1dHJKGI1SI4dNSAk+II8TWtvaoQ598OMq6F+Dfz8dPjVOCgta/+9eUwBloi0qppbjQEubP+rXc+jf9HX3sae58F+la3H1KskIiK5qjOTUcRKAyyItPYls7x0wK8/gP57wYHXwRY/yXSpskLKAiwzuw84BPjGObezv68S+COwwj/tYufc1FSVQUQSM3n6ojbB1ZLuxwU2jnkMBh6U1jKJiIikSsy0+HaCofbSAPO2AXL1MnjxUvjgSei5DRz1IOx4OJhlumRZI5U9WA8AtwIPhe2/yTl3fQrvKyLtiGdQbhca+V/3kwJv+tN7sNl2GSitiIhIanRmMoqCSQNs0VgPb90Kr98AzY2w14VeVkvXHpkuWdZJWYDlnHvdzPql6voi0jHxDMrdgpW82/3M1vfs0+1xZii4EhGRPNOZySgKIg2wxX+nwwvj4ftPYeAhMGoCbNIv06XKWpkYg3WWmZ0EzAbOd86tzEAZRApGeG/VD+sbYw7KffrpKTxUfEXrsZ80Pc7EA3ZJd7FFRERSrrO9UHmbBtjiu0+8wGrxi7DZDnDC07D9vpkuVdZLd4B1O3AV3hj6q4AbgN9FOtHMxgBjAPr27Zuu8onklUi9VdEsq6ml4oenqPCDq1rXlf1KpzAxX1vjRESk4BVUL1Qi6tfCG9fDW7dBcTfY/2oYdip06ZrpkuWEtAZYzrmvW16b2d3AczHOvQu4C2Do0KHh4+5FJA6RBu9GM6XHtfDS+97GoN9QeuTdzEph2URERLJB3vdCJcI5+OApbxKLNcvgp8d6MwdvtFWmS5ZT0hpgmVlv59xyf/MI4IN03l+k0MS7YvyS7sdBs/d67uArGVJxTgpLJSIikjs6uk5WzvlqAUwbB0tnQe/B8JsHYZthmS5VTkrZ8spm9hjwFjDAzL40s98D15nZAjObD+wNnJeq+4tI9EG6ZaUllJeVYjSHTMN+SP3VHDdnAFVzq9NVRBERkazVkmpfXVOLIzAxVF7Vkz98D8+fj7tjL2qWzufiht/zy+8voerbPAwi0ySVswgeG2H3vam6n4i0FW3wbuVhO1ExcAO4dtvW/T+tu4tVbAhxrv8hIiL5rWB6bmLozDpZWa+5Cd57EF65CldbwyPNI7lu/ZGsZkNYtT7uBZelrUzMIigiaRJ18G7v7+HaHVvP61/3MC6oQzve1EJQBSwiko9iLelRSM/4zqyTldU+fxumjoWv5sO2wzl5+a95ffWWIafkTSCZAQqwRPJcm8G78x6DO05r3exX92ib98Sz/geoAhYRyVd53XOTgM6sk5WVVi+Hly+H+VNg43L49X2w02jeuGhqxNNzPpDMEAVYIoWksmfgdd9fUDXkHko7sf6HKmARkfyUtz03CersOlmZ1pJlsqJmDeds+DJj3JOU0Ai/PN/76boBkIeBZIYpwBLJA3Gl6QUHV78aD3tfRIW/2dEUP1XAkuvMbBvgIWArvLk073LO/dXMNgWmAP2AJcBvnHMrM1VOkXTL9S/cyUpfD0+171laghmcN2Uek6cvyuq0+JYsk2FN7/FQ14fYrnE5M9xuNI6cwP577hFybq4HktlGAZZIjosrTS84uNp/AvzirNbNzqz/kesVsAjQCJzvnHvPzDYC5pjZS8ApwCvOuUlmNh4YD4zLYDlF0iqXv3AnO329pZ7MtbT4R6a9xv9xDyO7zuHT5q04Zf2FzGweTPmbtey/Z+i5WnA5uRRgiWSJjra2xUzT23kzmBA0aPXU16H3T5NW5lyugEUA/LUZl/uv15jZQqAcOBwY4Z/2IDATBVhSQHL5C3eq0tdzJi1+/Q/w5o08XH8zjUVFTGo4hvuaDmQ9JUD0LJN8XXA5E5NxKcASyQKdaRWL9qAsWrUUJhwR2HHJN9ClW3IK7MvlClgknJn1A4YA7wBb+sEXzrnlZrZFBosmkhG5+oU7VenrWZ8W7xx8+E948VJY/SWvFe/FpT/8hq/ZNOS0QsoyyVSvowIskSzQmVaxSGl6pxc/y7iSxwM7KlclrazhcrUCFglmZhsCTwHnOudWm1m87xsDjAHo27dv6gooInFLVfp6VqfFf/0RTLsQlrwBWw2CI+/mh++3ZfXTC6CAs0wy1etY1P4pIpJqsVrFquZWM3zSDPqPf57hk2a0WT1+7KgBlJYUt26/1HVs2oIrkXxgZiV4wdUjzrmn/d1fm1lv/3hv4JtI73XO3eWcG+qcG9qrV6/0FFhEYgqvFyE5gUWqrtsptTUwbRzcsSd8/QEcfAOMeQ22/QUVQ8qZOHoQ5WWlGFBeVsrE0YMKqlE0U72O6sESyQLRWsXKepS027UdnKY3qy4oJbCoBC77NsUlF8lt5nVV3QssdM7dGHToWeBkYJL/5zMZKJ6IdECq0tezKi2+uRnm/h1euQJ++B6G/hb2uRR6hKYDFnqWSaZ6Hc05l9IbJMPQoUPd7NmzM10MkZQJzxEGr1WsW5ciamob2pxfXlbKrPH7hO4MninwZ3+Eg69PVXFFUsLM5jjnhqb5nnsCbwAL8KZpB7gYbxzWE0Bf4HPgKOfc97GupbpKRNLiy9kw9QJYNhf67gEHXpvUCazySbTvVx3tyYu3nlIPlkgWiNYqdt6UeRHPD+nabmqAqzYPbB//FFVrf8LkSTMy38ImkuWcc28C0QZc7ZvOsohI6mRiJrmkW/sNvFwJ8x6BDbeC0XfDoKMgzjGjhShTvY4KsESyRKRu/MnTF8Xu2v7+M/i/wYEDFyymanFDTq3TISIi0iIVgVCurV/VRlMDvHsXzJwEDbUw/FzY6wLotlGmS5YTMpEmqUkuRLJYzAG18/8RGlxdthI23CLmjDkiIiLZqiUQqq6pxREIhMInd0pUTteLn7wKtw+H6RfDNrvDGW/DyCsUXGU59WCJpFkirXNRu7b/Ow4W/itwYtBMgVm/ToeIiEgEqZpSOyfrxZVL4cW/eHX9Jv3g2MfhxwcoHTBHKMASSaOOpCm06doOnswC2kzDntXrdIiIiESRqkAop+rFhlqY9Vd48yawItjnEtjjT1DSPdMlkwQoRVAkjTqdphAcXP1o74hrXEVLK9x7YK+Y62mJiIhkUrSAp7OBUFauXxXOOfjoWbh1GMycCAMOgrP+A3uNVXCVgxRgiaRRh1vnnAsNrg67BU6qinhqpIUFj9ytnKfmVCc9r11ERCRZUhUIZf2Cu998zDd/OxCeOJGPVzrOKrmKqu2vhp5bZ7pk0kFKERRJow6lKfzwPVzXP7B91mzYfIeY9wlPKxw+aUZK8tpFRESSJZVTamflgrt1q2DmtTS/cyfdmrtxeePJPNy0H03ri3kll2Y5lDYUYImk0dhRAyIueBe1de6z1+HBQ1s3nz1kHtfe8xnLav6bUMWTkwN8RUSk4IQHQlVzqxmeb+s6NjfD+496a1qt+5Zni0dyZe1ovmfj1lPUCJrbFGCJpFFCrXMvXeYNdPVVHf5Rh9fxyKkBviIiIuTB+lWRVL8HU8dC9WzY+mdw3BOcd8tyXIRT1Qiau6IGWGY2FTjDObckfcURyX9xpSlc2w9qVwa2K1cxuRNpfgn3nInkCNVVIpmTikWBg6Vq2vaMWPet12M192HYoBdU3AG7HA1FRfQpm5FTjaCp/nvPB7F6sB4AXjSzB4HrnHMN6SmSSO7r1MMneDKLjbeGP38IdC7NL5V57SIZ9gCqq0TSLh29S3mR3t7UCP+5B169BhrWwR5nwq/GQfdAOmAuNYLmZa9iCkQNsJxzT5jZ88BlwGwz+zvQHHT8xjSUTyTndOrhExxcDTkRDr+1dbOzaX5ZOcBXpJNUV4lkRjp6lzKV3p60HprPXodp4+Cbj7ylVQ68Fnq1DZpyqRE0r3oVU6i9MVgNwDqgG7ARQZWWiETWoYdPQy1M2Cqw/ZuHYMfDQ07JpRYukTRTXSWSZunoXcpEvZeUHpqaL+DFS+CjKijrC0c/DAMPAbOob8mVRtC86FVMg1hjsA4AbgSeBXZ1zv2QtlKJ5KCWFq9IrW0Q4+Gz9C24/4DA9oWfQY9N25yWSy1cIumiukokM9LRu5SJeq9TPTQNdfDvW+CNGwAHIy6G4WdDSXaOpeoITZoVn1g9WH8BjnLOfZiuwojkqvAWr0giPnymXgjv3hnYvrwmL1q4RNJIdZVIBqSrdynd9V6Hemicg0VT4YWLoGYp/OQwGDXB673KM8qmiU+sMVi/TGdBRHJZpBavYBEfPsHjrYD+dY/S59pX1SslkgDVVSKZka9ZFQn30Hy72Btn9ckr0GsgnPQM/GhESsuYSfn6955sWgdLJAlitWyVR3r4hAVX/eoeBTQbj4iI5I58zKqIu4emfg28dh28fbuXAjhqIgz7IxSXpLnE6deZv/dCmeJdAZZIEkRr8SovK2XW+H1CdwYFV9/Tk13rbg85rNl4REREMqPdHhrnYP4T8NJlsPYrGHIC7FsJG/bKXKFzRCFN8a4ASyQJ4mrxam6CK4MmrzjqAXb7e9eI16uuqaX/+OfzunVHRESkM1LVGxK1h2bZPJh2IXzxDvTZFY55BLYe2un7FYp0TfGeDb1kCrBEkqDdFq/vPoFbdg284Yx3YIuBUVdvB3Dkd+uOiIhIR6W1N2TddzDjKpjzAPTYDA67FQYfD0VFyb1PnkvHFO/Z0kumAEskSaK2eL1zF0wbG9i+9NvWHO1IPV/hlDIoIiISKi29IU2NMOd+mHG1N+Zq99NgxHgoLUvO9QtMOqZ4z5aFkBV6i6TSLUNDg6vKVSEDYCuGlHPkbuUUx5iaHbSAn4iISLCU94Ys/TfcNQKmXgC9d4HTZ8GBkxRcdcLYUQMoLSkO2ZfsKd6zZSFk9WCJpErYTIFUrmpzStXcap6aU02TczEvpQX8REREAlLWG7J6Gbx4KXzwJPTcBo56EHY8POYalRKfdEzxni0LISvAEkmFOIIraH/9LNACfiIiIuGSvuBtYz28dRu8fj00N8KvxsHwc6FrjySVODHZMFFDKqR6av9sWQhZAZZIkKQ80IKDq8HHQ8Xfop4aq8vaIK8eqiIiIsmS1N6Q/06HF8bD95/CwENg1ATYpF9yC5yAbJmoIRdly0LICrBEfJ1+oNWvgYlbB7ZPrILt9o75loTWzxIREZFWne4N+e4TeOEiWDwdNtsBTngatt83eQXsoEKazjwVsmEBbAVYIr5OPdCq58DdQQHRuKVxDYQdO2oAY598n4amwBiskmJTSqCIiKRFvn7Jjql+LbxxA7x1KxR3hZFXeTMEdom8NmW6FdJ05vlKAZaIr8MPtDduhFeuCGxHGW8VVfj8FrHnuxAREUmKgvuS7Rx88JQ3icWaZbDLMTDyCthoq0yXLEQhTWeerzRNu4gv2oMr5gPt1mGdCq4mT19EQ3NoRNXQ7Jg8fVFC1xEREUlUrC/ZeeerBfDAwfDU72HDXrz+y0cY/t+j6T9hDsMnzaBqbnWnb1E1t5rhk2bQf/zznbpmIU1nnq/UgyU5KRUpDQnPPBM8mUVxV7h0RcL31ANOREQypSDqoB++h1cnwOz7oHsZHHIzVUX7ctE/P6K2wfs9k9Fzl8zewEKazjxfKcCSnJOqlIaEHmjBwdWwMXDQ5A7dUw84ERHJlLyug5qb4L0H4ZWroK4GfvYH2PtiKN2EyZNmJD09Ltkpd4UynXm+UoAlOSeVecPtPtAa6+HqLQLbJzwF2+/X4fvpASciIpmSt3XQ5+/A1Avgq/mw7XA48DrYaufWw6noucu13sBsmc48XynAkpyTsYfYVwvgjj0D2+fMh022TegSkVIbJ44epAeciIikXS5+yY45RGDNV/DS5TD/cdi4HH59H+w0GsxCrpGKnrtc7A3MhunM85UCLMk5GXmIzbwWZl4T2L5sJRQlNkdMtNTGiaMHac0rERHJiFz6kh2tHrXm9Rxe9yy8dh00rYdfnu/9dN0g4nVS0XOXt72B0iEpm0XQzO4zs2/M7IOgfZua2Utmttj/c5NU3V/yVzpm1wkxcZvQ4KpyVcLBFRTYbE0iOSJKXVVpZtVmNs//OSiTZRQRT6R69GdN7zH4uYPhpcug355wxtuw72VRgyvwgsqJowdRXlaKAeVlpUwcPajT47iTfc14JGvmQkmuVPZgPQDcCjwUtG888IpzbpKZjfe3x6WwDJKHUp3SEJx+8Fn340IPJrrGVZBcy88WKRAP0LauArjJOXd9+osjItEE15fb2Ndc2uVh9i+ew2dNW8IJ/4Af7x/3tVLRc5fu3sCCW8csh6QswHLOvW5m/cJ2Hw6M8F8/CMxEAZZ0QKoeYsEPqyVJDK4gN/OzRfJdlLpKRLJQn7JSvqup4Ywuz3Bq8fM0UsSkhmN4YcPRzEwguMoX0TJjzn/ifc6bMi8nxtTlq3QvNLylc245gP/nFu2cL5IWLV3s506ZR21DY0hwdWfjwQzv/s9O3yPtqY0i0hlnmdl8P4Uwajq7mY0xs9lmNnvFisTXwhORODnHTTt/yoxuF3B2lyqmNf+Mfepv4MGiIzj3gJ3bf38eipYB0+QcjkCPltIG0y9rJ7kwszHAGIC+fftmuDSSz4J7rcpYw7zup7Yeq6i/knlueywJaXy5OFuTSIG6HbgKcP6fNwC/i3Sic+4u4C6AoUOHunQVUCRczNn1ct3XH8G0Cxm25A1qygZyeu2feWF1f/qUlTIxn37PBEXLjAkWvIxNXv8byTLpDrC+NrPezrnlZtYb+Cbaiaq0JF1autj3KPqQx7pOaN0/sO5+6ugGJC+NL5dmaxIpVM65r1tem9ndwHMZLI5Iu3JtLE7cX/Rra2DmRHj3bui+MRx8A2W7/Zbbi4rbnluAIs1cGMmymtqc+zeS69KdIvgscLL/+mTgmTTfXwpIvDPrLKupZVyXx0KCq351j7YGV0rjEyksfgNgiyOAD6KdK5INcmmW2pYv+tU1tdHT2Jqb4b2H4Jbd4J07YbeT4U/vwc/+AAquWoXPXFgctt5Xiz5lpTn1byQfpKwHy8wew5vQYnMz+xK4HJgEPGFmvwc+B45K1f2lsCXSUjO3+6mUsaZ1u1/do62vy9WFLpLXotRVI8xsMF6K4BLg1KgXEMkCuTRLbawv+hVDyuHL2TB1LCx7D7b5OZz4NPT+aYZKm/2CM2PCv/tAoJH4vCnzIr4/G/+N5INUziJ4bJRD+6bqniIt2n2At6jsSZn/crnblD3qbwW8B1I61q8QkcyKUlfdm/aCiHRCtLE4Dhg+aUZWNRRG+0K/vuYrqDoD5j0CG24FR9wFu/wGovTKSFuxxnpPnr5IMxmnUdZOciHSGXG15lX2bH15F0dyTf2RAGzSo4TLD90payojERGRWGKNxcm2sTbhwWAXGjm5+EXOK3ka5jfA8HNgr7HQbaMMljJ3RRvrHenfiIZApE66x2CJpEW0Fpk+ZaVQvzYkuDq+6QquqTuydbuuoTnl5RMREUmW4LE4kWTTWJvgJUuGFy1gWteLuLTkYdZtsSuc8RaMvFLBVQqEj9cqLytVpk4KqQdL8lK0lpqJQ9fCxMDD5NBu97FgVfeQ90ZMJRQREcliLT0X/cc/T6Spl7NlrE3FkHJ6/FBN11cuZUTz21Tblrw97DZ+fsDxSgdMMc1knD4KsCQvRcpDfqj3k2z35iOBky6v4YOLpkZ8f6yKSOtIiIhItoo2Hisrxto01MKsv7L/mzdBcRGMuITyPf5EeUn39t8rkkMUYEneCmmpqewJnwUdrFwFJF4RaR0JERHJZlk51sY5+Pg5mH4x1HwOO42G/a+CnltnrkwiKaQxWJL/gsZbedurWl8G54IHW1ffGHHdLK0jISIi2SzrxtqsWAR/r4ApJ0DXDeHk5+Co+xVcSV5TD5bkt+DgaoNeMPZ/IYdbKpwr/vUhK39oaN1fU9sQsWcql9YaERGRwpQVY23qVsNr18I7d0DXDeDA62Do76FYXz0l/6kHS/JTc3NocLVfZZvgqkXFkHJ6dG37wI/UMxVzdkIREZFC19wMcx+BW3aDt26DwcfBn96D3U9VcCUpUzW3muGTZtB//PMMnzQjYhZSOulfuqRVWiaIWL0cbhwY2D53AZT1jfmWeHumsjK3XUREJBtUvwfTLoQv/wPlQ+G4KVC+a7tv0+RR0hnZOD5eAZakTVr+AyyaBo8dE9i+9Lu4Wszinewi1irpIiIiBWntCnjlCpj7sJeOX3E77C9uqsgAACAASURBVHIMFLWfKJWNX47TQUFl8sQaH68AS/Jeyv8DPHs2vPdgYDtoMov2JNIzlRW57SIiIpnW1Aj/uQdevQYa1sEeZ8KvLoTuPdt/ry8bvxynWqEGlamSjePjNQZL0ial/wEqe4YEV/3rHk0oBzfrZl0SERHJZp+9AXf+El4Y56UBnv5vGDUhoeAKsvPLcappRuLkysbx8erBkrRJ2eKHQZNZfOj6c3D9BCDxFiH1TImIiLRj1Zfw4iXw4T+98c1HPwIDDwazDl0uqxdGTpFCDCpTKRvHx6sHS9Im0ppTnf4PEBRc/V+XU1qDqxZqERIRkUKXlBnWGurg9clw68+88c4jLoYz34WfHNLh4ApS9N0gy2Vjj0suy8YsJPVgSdokdYKI+jUwMWiRwhOruOnuHyKeGqlFSINLRUSkEHR6vI9zXkA1/SJYuQR+cpiXCtjO7LzxKsTJo7KxxyXXZVsWkgIsSauk/AeongN37xPYHrcUSsvoUzYjrjQDDS4VEZFC0alJJL5dDC+Mh/+9DJsPgJOegR+NSHoZs+3LcaoVYlBZaBRgSdaK2Mu0doo3FWyLoJkC420RKsQZi0REpDB1aLxP/Rp47Tp4+3YoKYVR18CwMVBckqJSFp5CCyoLjQIsyUqRepl2rtof7MvASWHTsMfbIqTBpSIiUigSmkTCOZj/BLx0Gaz9CgafAPtdDhtukYaSioTK5eEcCrAkK4X3Mi3pflzgYFEJXPZtxPfF0yJUiDMWiYhIYYp7vM+yeTDtQvjiHT6y7flL/ZV88/Egxm7bQMWQNBdaCl6uD+dQgCVZKbg3KTi4eqhxJCdd/WTrdkdaNzS4VERECkWk7I69B/Zi8vRFnDdlHgN7NvC33s/Tf8k/qO+6CVc2n8aj6/fEUQQ59qVW8keuD+dQgCVZqU9ZKV/XrOF/3U9q3Xfy+nH8b+Of07Kno60bGlwqIiKFJDi7o6XurG9o4PjiV7ig7gk2/KyWT7Y7gTOqR7FofegKPrn0pVbyR64P51CAJVnp8j17sP/LR7Ru/6zub6wt2YyJQb1MnWnd0OBSEREpRJOnL2Lnxg+5ouuD7Fi0lH837Uhl48msW/Zjlq1K3pfaXB4/I5mX68M5FGBJ2rX70J3/D/Z/+Q+tmz+qe5jeZRswMey8XG/dEBERSavVy7hw3WQO7/Zvqt1mnLH+bKY27w4Y5tfJyfhSGynDZOw/3ueKf31IzQ8NCrikXbk+nEMBlqRVu2l9U06Ehc8G3lC5ik+jXCvXWzdERETSorEe3roNXr+eA4rX89fGI7ij8VBq6d56SkvQE+tLbby9UpEyTBqaHSt/aAByb8ICSb9cH86hAEvSKmZa3zM7hp4cNg17uFxv3RAREUm5/073Fgv+/lMYcDCvbXM2d0yvoZa2dWesL7WJjHuOJ5NEY7ukPbk8nEMBlqRVtIfurLrAeCt+tDecVNVuS1mut26IiIikzHefwAsXweLpsNkOcMJTsP1+7A9M3CB6/RrtS20i456jZZiEU0q/5CsFWJJWbR+6jiXdjw9sHvp/sNvJcbeU5XLrhoiISLzinjRi/Tp4/Xp461Yo7gojr4LdT4MuXVtP6Ujdmci450gZJpEopV/yVVH7p4gkz9hRAygtKQagjDUhwdW+9ZMZ/tI2rZVItJYyERGRQtLS6FhdU4sj0OhYNbc6cJJzsOBJuGUovHkj7Hwk/GkODD87JLjqqGjBUKT9FUPKmTh6EOVlpRhQVlpCSbGFnKOUfsln6sGStGppMXt56lPc2nBp6/4f1z3IekpaFzWM1uqldAIRESk07abnffUBTLsQls6C3j+Fox6AvrsntQyJjnsO7yXTtO1SSBRgSdpVdH+PiqDgql/doyHHaxuaKDajybk271U6gYiIFJpojYvralbA8xfA7HuhexkccjPsehIUFSe9DJ0d96yUfikkCrAkvb5aAFNOaN0MD65aNDlHaUmxZggUEZGCFz5+uYhmji5+lbFdptD07jr+2eVASkdcysFDd4xxlc5TkCQSH43BkvRZ/BLcd4D3evTdVB3+ERbl1PKy0pD87ZZtPdjl/9u78/ioqvOP458nk4WASFAEJQquFatU0Ui1WMUVtbbGrUq11S5SW21dEEXlB1gXqGitVm2llrrvyKIgiEWq4gYIdQVRRAkg4BIUSMhk5vz+mEkyk8wkkziZO5n7fb9evpg7uTNz7p14nzz3nPMcERG/iZ2/fKB9wLTCUYwr+BfL3c6cVHMjl28+h8uf+TR+TpaIeEY9WJIZCydFhjH02hd+9hhs25sJ4+fSdBAgGNQPO1BCJSIiflc+oJSi6g0wZwwnhOex1m3HH2ou4unwoRC9Val1pUSyhxIsaV/hMDw/Bl65nVfyDmLYyt/T7a6ljBjiko4pd2hldxER8Z+EhSD67wCv/50T/nsTWA0cdhnHPL8fm+nU5PUqBCWSHTREUNpPsAqePA9euZ1Hwsfy8y2XsIni+vKyJZ0LEr6sVIUsRCRDzGySma03s3dintvOzOaY2fLov929bKP4Q6JS7NOevJ+Pr98f5ozm5dp+zBk8DY4ZQ0lJ4l9JFYISyQ7qwZL2sflzeOQsqFjI3/LP45ZNx0LMjKuqYIii/DwVshARr90L3AHcH/PcSOA/zrnxZjYyun2lB20TH6jrtYotYrGLreP/8h/kuMAiVoR35LzaEcwLD6D4ua8Zt83qVpdM9wuVgpdsoQRL0u/z5fDQ6fDNZ7xx8K3c8lLPhLttrApy65kH6GIoIp5xzr1oZrs2evpkYHD08X3APJRgSTuo67WqS5Q6sZXf5U/ngsAz1JLH+OBZTAqdEFknkoZ5VvNHHgW0vWR6Lmp8LutGy4CmHUjmKcGS9Fo5Hx79GeTl899D/80F8/KAxIsG9y4pViELEclGvZxzawGcc2vNLPFdIsDMhgHDAPr06ZOh5kmuaFhA2HFi3utcU/AQpfYFU0M/YFzwZ6xjuyavqZtnpfgZr8XFmEUySAmWpM9bTxCa+jsqXE9+Xn05q5+HkEucXGkog4jkAufcRGAiQFlZWaLCqCJJramsYi+rYGz+fQwKvMt74b5cUnMhC1y/pK/RPKvEkhX4UOEP8YISLPn2nIMXb4YXrmeR+y7nb72EjWwDCYuwR2hNKxHJYuvMbKdo79VOwHqvGyQ5qKqSP3d5mFNrZ7KJYkYFf8kjoaMIEVnvqqS4gM01tQRDDbFUNyeTa7wYc+zzIpmmKoLy7VRVwrUl8ML1zMo7gnO2XhlNrpIrjQ4NFBHJUtOBc6OPzwWmedgWyTXhMLx5P/ztIM4IzeBJdzSDt/6FB0PHEiJAcUGAv555AEvGHMeE0/entKQYIxI7dXMyudjFmOsoIRWvqAdL2l5157O34R+H1W/+bsswXEylwER0sRORbGJmjxApaNHDzCqAMcB44HEz+zXwKXCGdy2UnFKxEGZeDmsWwy6HYD9/ik6f9aDL7GVsTBCDNc8qdXXnSYU/JBsowfK5Nlfd+d+jMOW3DdtjN9J7/NyE3fMBM8LO6WInIlnHOTc0yY+OzmhDJLdtWg/Pj4UlD8E2O8Kp/4T+Z4AZ5Tupyl26KCGVbKEEy+faVHVn2kWw+IHI410OgV/PBki6LoeGNIiIiC+FgvDGRJg3HoJVMOhiOHwEFHX1umUi0o6UYPlcsuo6qyurGDR+btNu9pv2gC2fR3Y6chQcMaL+NeqeFxERiVoxD569EjYshT2PgePHQ4+9vG6VeEgLIfuHEiyfS1Z1x6D++bphg+XTvtuww8+nwB5HNXmduudFRMTXKj+F2dfA+9Oh+65w1iOw9wlgzc9RltymhZD9xZMqgma20szeNrMlZrbQizZIRKKqO0Z8gfUiang/cGbDE5e8kzC5EhER8a1gVWQo4B0Hw/I5kVEev38d+p2o5EqanZIhucfLHqwjnXOfe/j5vtFcl3SiYX2xPVo72wZeLrq44c1GrYf8ooy2X0REJGs5B0ufgdlXR3qv9j0Fjrseuu3sdcski2ghZH/REMEcl0qXdONhfYOi1QAH5y3m3sIJDc93msJ8JVciIiIRG5bBs1dE5lv1/C6c+zTsdrjXrZIspIWQ/cWrBMsBz5mZA+52zk30qB05ry1VAkcM2Zv9ph7LnrYagM9cd44M/4NxCdav0oRNERHxm2feWMY3s6/n9NoZVFknVva/mu+VD4eA7ltLYskqLWtt0Nzk1ZVgkHNujZn1BOaY2VLn3IuxO5jZMGAYQJ8+fbxoY05oS5d0+bTvUrde8NeuM6cVT2JcgsRJEzZFRMRXwmHefPouvv/mTWzP1zwWGsyE2jOpWtKdcbuvU+yTpFRp2V88SbCcc2ui/643synAQODFRvtMBCYClJWVuSZvIilpdZf02G4Njw/+Ddv+6BbmJ3nvNq2hJSIi0hGtfhOevYIDKxaw2O3Jr4IjeNvtHvmZYp+kQJWW/SPjCZaZdQHynHPfRB8fB/wp0+3wi5S7pENBuK5Hw/bZT8JexzZ5v9ghgcmyXk3YFBGRjq4u3lVXrmNMlyf5ceg/WJcduDz4WyaHfohrVIhZsU9E6njRg9ULmGKRkqX5wMPOuVketMMXUuqS3rQBbt6zYXv4B9C1V5P3ajwkMBlN2BQRkY4g2TziqYtXM+qpJZwWns1lRU/QuXYr97oT2eGI0bw6dw1OxQpEpBkZT7CccyuA/TP9uX7WbJd0xSK4J2ZNq9FfQV7i5dESDQlsTBM2RUSkI2huHvHzMyfzpE2kX8EqXgrtx9jac/nIlRKYsoKQc03Wi2xN7FNxKJHcp3I3OaJNF+wF/4IZlzVsj93Y7O7NDX8wUKAQEZEOI9FNw5Lgero+fT53hOdTYT34bc2lzA6XUVf5KeQiaZWLPuOA0gSxr7meMRWHEsl9SrByQJsu2I+fC+9NjTzeawic/XiLn5OsYEZpSTHzRx6V4BUiIiLZKfamYRE1/CYwkwvzp5EXCnNPwVlM2HQ8WylM+vq65Kpx/GsuJqs4lIg/JB4LJh1KcxfshMZ2a0iuhoxLKbmCSMGM4oJA3HMaEigiIh1RZM6U45i8RTxXeAUjCh5nXnh/zi66gx4/Gk1eQctzqhKN7GguJrdl6RQR6XjUg5UDUr5gOwfXltRv/q7wBoZ0OpnyFD9HaziIiEiuGDuoiKLnr+VwW8LycCln11zFG3yPbWrzufSxJXQrLqBTQR6VW4LkmdUPD4xVV9gi1Qq7rV46RUQ6JCVYOSClC/aXH8PtB9RvHr11Ah9VlzKvlWO/tYaDiIhksxbnJG/9Bl6cwLGv3kWwoIjb+BV3bBpMl+JiqKnlqy1BACqrghQXBLj1zEjsTLbkSWsq7Ka8dIqIdGhKsDqo2ADSrbiAgoARDDXcN4u7YL/9JEz+df3P+lX/m2qKAI39FhGRjq1xPNxcU1sfD+PmJB/QG95+Ap77P9j0GRxwDgXHjOHibXpyMTBo/Fwqq4Jx710XI+vmWSVK3AaNn9ticmXRtkyYvYzTDirlhaUbNBJEJIcpweqAGt8tq6wKUpBndO9cQOWWYPwF+74fw8cv1r921+qHm7yfxn6LiEhHlCgeNlYVDDH12Wcpf/MJWPUa79meXLP1T6xf2p8RfYOUD4js19Jw+2QjOFqKobEl3VdXVjF50WrGndpfSZVIDlOC1QElmkAbDDs6F+azePRxDU+O7Ra3z6BOU6BaY79FRCQ3tLQ+YwnfcHn+4wzdOpet67pzXfi3PFTzQxx50KjiblvnRzVXYRdo8jONHBHJfaoi2IFMXbyaQePnJryQQ6O7aI2SK8ZuVBVAERHJKcl6j/IIc05gDvOKLuOswAtMzj+Rk+02Hqw5IpJcRcVW3G1rjGzudaoaKOJPSrA6iLphEMmSK4i5yxabXO1xdP0CwuUDShl3an9KS4oxInfXNExBREQ6qkS9SwfbUp4pvIbrC/7Ne+G+nBr+M4UnTWDZxkCCd4gfAtiWGNnc65L1fmnkiEhu0xDBDqKlYRDFBQFGHr1LfHJ1xr2w7ylx+6kKoIiIdDTJKgPGVuXrxZdcXfAwJwdeYS3bc2HNH1nSdTAjju9H+YBSJsxe1uIQwLbGyGSvU9VAEX9SgtVBNDecoLSkmD8dGuDoGWUNTw7/ALr2qt9ssWytiIhIFmpcyGJ1o7lTeaEaVs+6mV8EnyDfwiz9zu/od/po7izsHPc+XiQ7Wj9SxJ+UYHUQzU2inX/IAnjhxoYnx1SCWf1mS8FJREQkWyUawVFfKKLLu/zklZFQ+xH0+xEMuYF+2+2W8H28SnY0ckTEf5RgdRDJ7rzNrz4F5sXsGJ1vFavZ4KSLvoiIZLFEIzj62meM3vwAPLwYtt8TzpkMex6T8PUawZEZOs8iDZRgeSzVC1KiO2/zq+PnVyVKrqDltT1ERESyVewIjs5U8/v8aZwfmEGQfMYFhzJr0ylc+s0+lCd4rUZwZIbOs0g8JVgeau0FKW6YQYIy7Mm0dW0PERERr0VGcLzFMaGXubrgYXayL5kcOow/B4eynu6wsTZp7NQIjszQeRaJpzLtHmrugpRUqDYuufq47xkM6jSF3UbOYND4uUxdvLrJS7T+lYiIdDR1az/e/fh0Hghcy98K7+ALty0/rRnL8ODvI8lVVLLYqREcmaHzLBJPCZaHWn1B2rAMrtu+fvOlQfdz4orTWV1ZhaOhB6xxkqX1r0REpCOZung14596lWGb7uKZwqvY3a1idPh8Pix/mgXh7yR8TaLYqXWoMkPnWSSehgh6qFVD9168GeZe17B9zTpG3vIKVcH41yfrklcVIxERyWZ1c5I/q9zMWYF5zMx/lG5s5sHQMfyl9gw2sg3/ee7DVsVOrUOVGTrPIvGUYLWDVAtXpHxBuqE3BDc3bEfnW6lLXkREckHdnOR9at/nH4X30j9vJa+H+zE2eC7vu771+62prOLWMw9I+Y95rUOVGTrPIvGUYKVZawpXpHRBaqaYhYpXiIhILpg061Wu59+cVvQya912/KHmIp4OHwpY3H69S4pb/ce8RnBkhs6zSAMlWGnW2ko6zV6QWqgUqC55EZH2Y2YrgW+AEFDrnCvztkU5qLYGXv8HD1XfSGFekDtrf8KdteVsoVOTXWPjm/6YF5FspgQrzdI2bC8muVrX8zBO/Xo4a0bOiLtTpy55EZF2d6Rz7nOvG5GTPnwenh0JXyznrcBBXF11Np+4HeN2CZgRdk7xTUQ6FCVYafath+1t/gIm7F6/+drAv/HLV3vWF7NoPORQd/FERKRD+fJjmH0NLJvBKtuJMTUjeLNoIJvzaiHk6ncrLgio4q2IdEhKsNLsWw3be28aPP6Lhu0rPmb47UtSrhQoIiJp5YDnzMwBdzvnJnrdoGwUW9ipW3EBZlC5Jdi016lmC7x8K8y/jVryuC08lLtrjqeGAqgKUpBndO9ckPi1IiIdiBKsNGvzsL17T4KVLzVsp1gpMNWKhSIi0mqDnHNrzKwnMMfMljrnXozdwcyGAcMA+vTp40UbPdW4sFNlVbD+Z/UjLpyjvGghPDcKNq6C/U7n9OVDWLK5S9x7BcOOzoX5LB59XEaPQUQk3ZRgtYNWD9trY6XA1lQsFBGR1nHOrYn+u97MpgADgRcb7TMRmAhQVlbmmrxJjktU2CnWzrWfsPMz10P4bejVn5f638DIhV1ZvVHLjIhI7srzugG+1yi5mnrye3HbI4bsTXFBIO65uiGHzVUsFBGRtjOzLmbWte4xcBzwjretyj7JEqJt2czo/Pt5tnAke4ZWwIk3M23gQwz7b1HCm4Z1tMyIiOQCJVheikmu3gv3Zdfqh7nqqbeZunh1/fPlA0oZd2p/SkuKMaC0pLh+0q8WGhYRaTe9gJfN7H/AG8AM59wsj9uUdRonREaYMwLz+E/RcM4LzObR0JEM7XQXDDyfm+Z81Gxvl5YZEZFcoSGCXgjVwnXb12+etnUMi1wkqCQqYJFsyKEWGhYRaR/OuRXA/l63I1PaOp83trDT/vYh1xbcywF5K1gY/g7nBUeyIn8Pxh3fH2j+5l+p5hCLSA5RgpVpX6+Bv+xTv3lA9d1U0jVul1R7oNpSsVBFMUREJNa3mc9bPqCUourPCT0/lpNCc1lPd66xP/BwzSH0LunMaf12YMLsZVz62BLyzAi5ptPUSkuKmT/yqPQfmIiIR5RgZdKHz8ODp9VvHlY0mcrqrU12S7UHqrUVC1UUQ0REGmtuPm+zsSEUhDf+yQnzxoGrgh/8kZ5HXMENRV25gaYxJ1FypWGBIpKLlGBlypzRMP+2hu2xG7m8UfCB1geb1lQsbHMQFRGRnNWm+bwr5sGzV8KGpbDH0XDCn6HHXnG7JKswGDAj7JxGUYhIzsr5BCsrhsTdul9k7Q+Aom3hqsjjNq+Z1UYqiiEiIo0lm8+bZ8ZuI2fEx6bKT1n92HBK1z7HJ+Ge3FV0FYfuew7lPXZu8vpksSXsHB+P/1Haj0NEJFvkdIKVFUPiYsuwH/wb+NEtcT9u9ZpZ34KKYoiISGOJ5vNCw5C+1ZVVjH1qEXsvvZPvLL+H7UKOCbU/5Z7QiWytKWT6lHfArEksU8wREb/K6TLtnq4T5Vx8cnXmg02Sq0xrbk0tERHxp8bLgQTMYn7qGJK3gKftMvZZdifPBgdw1NZbuDNUzlYKgeRxVTFHRPwqp3uwPBsSF6yGG3o1bF/8FnTvm7a3b+uwx0wPSRQRkY4hdjTFbiNnALCHrWZM/v0cHnibpeFdGFpzDa+G9034+kRxVTFHRPwqpxMsT4YnfL0W/tKvYXvUBsgvTNvbf9thj5kckigiIh3Pnt0cZ2x+mF8GZlFFEWODv+CB0LGECCR9TbK4qpgjkvuyot5BlsnpBKst60R9KxWL4J7oWh69D4RhL6T9I1QJUERE2qLFP4LCYXjrMZ5mFIWBL3g8dAQTas/kC7olf1M07E/Ez7Ki3kEWyukEK6PDE5Y8AlMviDw+5EI4/sb0fwaqBCgiIq3X4h9BaxbDzCug4g06lZYxb487+dsbhXxZWUUgyQLBEFkkWHerRfxLN/4Ty+kECzI0POHZK+H1f0QenzIR9j+z3T5KVZlERKS1kv0R9M9Zb1C+6gV4837o0gNOvgv2H8rgvDzmRwdkNE7OINJrNe7U/r7+A0pEdOM/mZxPsNrd3YfD2v9FHp8/F0oPatePy/iwRxERyVqpzn1o/MdOgBDnBJ7nsuonYEkNHHohHHEFdGo6HFDFKkQkGd34T0wJVluFauG67Ru2hy+Drju2+8cq0ImICLRu7kPsH0GH5L3H2Pz76Je3igV5+3PwBROhZz+ao2IVIpKIbvwnZi7JuOpsUlZW5hYuXOh1MxpUfQV/3rVhe9R6yC/yrDkiIrnAzBY558q8bkdbZTpWDRo/N+Gd44AZYefibsBNXbya2556geE8wEmB16hwPbix9he8WngIlVW1ulknIm3mpyqCqcYp9WC11vqlcNf3I4+77gSXvQ9xizKKiIi0v2RzHOoKUtT1aOWFtlK+ZTInFd5CKBTi1uBpPFpwCl8SIFhVG7cv+Lvyl4i0nnq4m8rzugEdyrJnG5Kr/mfA8KVKrkRExBMtz3Fw/CD0Bt+bfjzMvZ71PQ+j6OKFXHrDJPI7dSEYih/BUlf5S0REvh0lWKl6cQI8clbk8fHj4bR7vG2PiIj42oghe1NckHjx391sLfcW3MS/Cm+hlgBn11zF0RW/YerKyMAVVf4SEWk/ngwRNLPjgduAAHCPc268F+1I2cNnwgezIo9/MQ12H+xla0RERJoUPcozo5Pbwh/yp/KrwEyqKeS64DncFzqOWvIh3LA2jSp/iYi0n4wnWGYWAO4EjgUqgAVmNt05916m29Ii5+CGHaG2OrL9x8Ww3e7etklERCSqfu6Dcyx8+m76LBpPT/uKJ2oP56bas9hASdz+dT1UqvwlItJ+vOjBGgh86JxbAWBmjwInA9mVYNVsgRt3ati+qgKKunrXHhER8aUWK3StfQtmjqBs1Wt81X0/hm25kueq+yR8r7oeKi35ISLSfrxIsEqBVTHbFcD3G+9kZsOAYQB9+iQOFO2m8lP4a/+G7dFfQZ6mq4mISGY1u9bV3sUw9zpYdC8Ud4cf3073AT9nYl5ek9dB0x4qVf4SEWkfXiRYicruNVmMyzk3EZgIkbVF2rtR9T55Bf59QuTxbofDuU9n7KNFRERiTZi9LC5JAtgaDLJ8xm3w3ONQ/TUMHAaDr4LihuGA6qESEfGOFwlWBbBLzPbOwBoP2tHUwknwzKWRx4ePgKNGedseERHxtcZV/cpsKdcW3Me+tZ/Azj+EE26CXt9N+Fr1UImIeMOLcW8LgL3MbDczKwTOAqZ70I540y5qSK5+er+SKxER8VzdnKmefMWtBXfyZNGfKLFNjCq4PDLCIklyJSIi3sl4D5ZzrtbMLgJmEynTPsk5926m2xHn9gPhy48ijy94GXbs3/z+IiIiaZSskMWVx+zG8uk3cYFNJp8wt9eWc6+dwugTD9ZC9yIiWcqTdbCcczOBmV58dpxQEK7r0bA9YgV02d679oiIiO8kK2TRa91/+ckHN0PeR7yUdzCjqn5GbbddGa25VCIiWc2TBCsrbP4cJuzRsP1/n0OgwLv2iIiILzUuZNHXPmMUD3Loa2/C9nvC2ZP54V7H8F8P2ygiIqnzZ4K19i24+4eRx9vvBX9Y6G17RETEt+oKWRRTzYX50zg/MIMg+YwLDmXWplO4dNM+lHvcRhERSZ3/Eqx3p8IT50YeH3gu/OR2b9sjIiK+1rukmOrKz3i66Bp625c8FTqM8cGhrKc7bKxtWPdKwwJFRDoEfyVY/7kOXro58vjHt8FB53naHBERkRFD9uaqp2qYHvoBc0IHscjtHffzqmCICbOX1SdYyQpiiIhIdvBPgnXvSbDypcjjX86Cvod62x4RERFiFwX+NasbrXtVp24YYbKCGLHvIyIiR4awrAAAB9dJREFU3vJiHazMe+ychuTqkneUXImISFYpH1DK/JFHURpd96qxuvWwGhfEgIYeLhERyQ7+SLB2OQSKtoWr10LJLl63RkREJKERQ/amuCAQ91xxQYARQyLDBte00MMlIiLe80eC9YOL4KpVUNjZ65aIiIgkVT6glHGn9qe0pBgDSkuKGXdq//rhf71b6OESERHv+WcOloiISCuY2fHAbUAAuMc5Nz4Tn1s+oDTpfKpIQYy344YJxvZwiYiI95RgiYiINGJmAeBO4FigAlhgZtOdc+952a6GghiqIigikq2UYImIiDQ1EPjQObcCwMweBU4GPE2woPkeLhER8Z4/5mCJiIi0TimwKma7IvpcHDMbZmYLzWzhhg0bMtY4ERHJXkqwREREmrIEz7kmTzg30TlX5pwr22GHHTLQLBERyXZKsERERJqqAGLX9dgZWONRW0REpANRgiUiItLUAmAvM9vNzAqBs4DpHrdJREQ6ABW5EBERacQ5V2tmFwGziZRpn+Sce9fjZomISAegBEtERCQB59xMYKbX7RARkY5FQwRFRERERETSRAmWiIiIiIhImijBEhERERERSRMlWCIiIiIiImlizjVZNzHrmNkG4BMPProH8LkHn5st/Hz8fj528Pfx+/nYwdvj7+uc67Cr9cbEKr//DqVC56h5Oj8t0zlqns5Py9pyjlKKUx0iwfKKmS10zpV53Q6v+Pn4/Xzs4O/j9/Oxg44/HXQOW6Zz1Dydn5bpHDVP56dl7XmONERQREREREQkTZRgiYiIiIiIpIkSrOZN9LoBHvPz8fv52MHfx+/nYwcdfzroHLZM56h5Oj8t0zlqns5Py9rtHGkOloiIiIiISJqoB0tERERERCRNlGCJiIiIiIikiRKsJMzseDNbZmYfmtlIr9uTSWa20szeNrMlZrbQ6/a0NzObZGbrzeydmOe2M7M5ZrY8+m93L9vYnpIc/1gzWx39HVhiZid62cb2Yma7mNkLZva+mb1rZhdHn8/577+ZY/fFd99e/Bw7UuG3+JIKv8eglvg5RqXKz7EsFV7EO83BSsDMAsAHwLFABbAAGOqce8/ThmWIma0EypxzvligzswOBzYB9zvn9os+dxPwpXNufPSPpO7OuSu9bGd7SXL8Y4FNzrmbvWxbezOznYCdnHNvmllXYBFQDpxHjn//zRz7T/HBd98e/B47UuG3+JIKv8eglvg5RqXKz7EsFV7EO/VgJTYQ+NA5t8I5VwM8CpzscZuknTjnXgS+bPT0ycB90cf3EfkfMSclOX5fcM6tdc69GX38DfA+UIoPvv9mjl3aTrFDWs3vMaglfo5RqfJzLEuFF/FOCVZipcCqmO0K/PWHhwOeM7NFZjbM68Z4pJdzbi1E/scEenrcHi9cZGZvRYdn5PywAjPbFRgAvI7Pvv9Gxw4+++7TyO+xIxWKL6nx1TWojXSdSsDPsSwVmYp3SrASswTP+Wks5SDn3IHACcCF0e558Ze/A3sABwBrgVu8bU77MrNtgMnAJc65r71uTyYlOHZfffdp5vfYkQrFF0kHXacS8HMsS0Um450SrMQqgF1itncG1njUloxzzq2J/rsemEJk2IvfrIuO2a0bu7ve4/ZklHNunXMu5JwLA/8kh38HzKyAyAX3IefcU9GnffH9Jzp2P3337cDXsSMVii8p88U1qK10nWrKz7EsFZmOd0qwElsA7GVmu5lZIXAWMN3jNmWEmXWJTgDEzLoAxwHvNP+qnDQdODf6+Fxgmodtybi6C3LUKeTo74CZGfAv4H3n3F9ifpTz33+yY/fLd99OfBs7UqH40io5fw36NnSdiufnWJYKL+KdqggmES3V+FcgAExyzt3gcZMywsx2J3JXESAfeDjXj93MHgEGAz2AdcAYYCrwONAH+BQ4wzmXk5Nskxz/YCJd5g5YCfy2bhx3LjGzw4CXgLeBcPTpq4mMzc7p77+ZYx+KD7779uLX2JEKP8aXVPg9BrXEzzEqVX6OZanwIt4pwRIREREREUkTDREUERERERFJEyVYIiIiIiIiaaIES0REREREJE2UYImIiIiIiKSJEiwREREREZE0UYIl4iEz28XMPjaz7aLb3aPbfc1slplVmtkzXrdTRET8qZk4dYSZvWpm75rZW2Z2ptdtFckWKtMu4jEzuwLY0zk3zMzuBlY658aZ2dFAZyLrMpzkbStFRMSvEsUpYDLgnHPLzaw3sAjYxzlX6WFTRbKCEiwRj5lZAZHANAk4HxjgnKuJ/mwwcLkSLBER8UpzcSpmn/8BpzvnlnvQRJGsku91A0T8zjkXNLMRwCzguMZBS0RExEstxSkzGwgUAh950T6RbKM5WCLZ4QRgLbCf1w0RERFJIGGcMrOdgAeAXzrnwl40TCTbKMES8ZiZHQAcCxwCXBoNViIiIlkhWZwys22BGcAo59xrHjZRJKsowRLxkJkZ8HfgEufcp8AE4GZvWyUiIhKRLE6ZWSEwBbjfOfeEl20UyTZKsES8dT7wqXNuTnT7LqBftPztS8ATwNFmVmFmQzxrpYiI+FXCOAVcBRwOnGdmS6L/HeBVI0WyiaoIioiIiIiIpIl6sERERERERNJECZaIiIiIiEiaKMESERERERFJEyVYIiIiIiIiaaIES0REREREJE2UYImIiIiIiKSJEiwREREREZE0+X+Q1ibgNDjlZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x396 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig, (ax1,ax2) = plt.subplots(1,2 , figsize = (12,5.5) )\n",
    "\n",
    "ax1.plot(x1, y, 'o')\n",
    "ax1.set_xlabel('X1')\n",
    "ax1.set_ylabel('Y')\n",
    "ax1.set_title('Contant variance of residuals')\n",
    "ax1.plot(x1, lr.predict(x1.reshape(-1,1)))\n",
    "\n",
    "ax2.plot(x2, y2, 'o')\n",
    "ax2.set_xlabel('X2')\n",
    "ax2.set_ylabel('Y')\n",
    "ax2.set_title('Non contant variance of residuals')\n",
    "ax2.plot(x2, lr2.predict(x2.reshape(-1,1)))\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quantile regression**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.DataFrame(data={'X': x2, 'Y':y2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'smf' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-17c9bc5331b1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mmod\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msmf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mquantreg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Y ~ X'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mres\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmod\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mq\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m.5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mquantiles\u001b[0m  \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0.05\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0.95\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'smf' is not defined"
     ]
    }
   ],
   "source": [
    "mod = smf.quantreg('Y ~ X', data)\n",
    "res = mod.fit(q=.5)\n",
    "\n",
    "quantiles  = np.array((0.05, 0.95))\n",
    "def fit_model(q):\n",
    "    res = mod.fit(q=q)\n",
    "    return [q, res.params['Intercept'], res.params['X']] + \\\n",
    "            res.conf_int().loc['X'].tolist()\n",
    "    \n",
    "models = [fit_model(x) for x in quantiles]\n",
    "models = pd.DataFrame(models, columns=['q', 'a', 'b','lb','ub'])\n",
    "\n",
    "ols = smf.ols('Y ~ X', data).fit()\n",
    "ols_ci = ols.conf_int().loc['X'].tolist()\n",
    "ols = dict(a = ols.params['Intercept'],\n",
    "           b = ols.params['X'],\n",
    "           lb = ols_ci[0],\n",
    "           ub = ols_ci[1])\n",
    "\n",
    "print(models)\n",
    "print(ols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'models' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-13-8a899a67c6ec>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m6\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m     \u001b[0myn\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_y\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmodels\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m     \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0myn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlinestyle\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'dotted'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'grey'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'models' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAFpCAYAAAC8iwByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEBRJREFUeJzt3V+I5fdZx/HP06yxWGsrZoWSP03ErXUpQusQK4JWWiXJxeamlgSKVkIXqqmgpRBR2hKvrIggROuqpbbQxuiFLrISQSOV0pRsqYYmJbCmtVkiZK01N6VNo48XM63jZHbnt5szu8/ueb1g4PzO+c6Zh+8O887vzJlfqrsDAMz1kks9AABwbmINAMOJNQAMJ9YAMJxYA8BwYg0Aw+0Z66r6cFU9U1WfP8vjVVW/X1WnqurRqnrD6scEgPW15Mz6I0luOcfjtyY5tPVxNMkfvvixAIBv2TPW3f3JJP95jiW3J/lob3o4ySur6lWrGhAA1t0qfmd9bZKnth2f3roPAFiBAyt4jtrlvl2vYVpVR7P5Unle9rKX/ehrX/vaFXx5AJjvs5/97H9098EL+dxVxPp0kuu3HV+X5OndFnb3sSTHkmRjY6NPnjy5gi8PAPNV1b9d6Oeu4mXw40l+futd4W9M8mx3//sKnhcAyIIz66r6RJI3Jbmmqk4neX+S70iS7v5QkhNJbktyKsnXkvzifg0LAOtoz1h39517PN5JfnllEwEA/48rmAHAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAw3KJYV9UtVfVEVZ2qqnt2efyGqnqoqj5XVY9W1W2rHxUA1tOesa6qq5Lcl+TWJIeT3FlVh3cs+80kD3T365PckeQPVj0oAKyrJWfWNyc51d1PdvdzSe5PcvuONZ3ke7ZuvyLJ06sbEQDW24EFa65N8tS249NJfmzHmg8k+buqeneSlyV5y0qmAwAWnVnXLvf1juM7k3yku69LcluSj1XVC567qo5W1cmqOnnmzJnznxYA1tCSWJ9Ocv224+vywpe570ryQJJ096eTvDTJNTufqLuPdfdGd28cPHjwwiYGgDWzJNaPJDlUVTdV1dXZfAPZ8R1rvpzkzUlSVT+czVg7dQaAFdgz1t39fJK7kzyY5AvZfNf3Y1V1b1Ud2Vr2niTvrKp/SfKJJO/o7p0vlQMAF2DJG8zS3SeSnNhx3/u23X48yU+sdjQAIHEFMwAYT6wBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGG5RrKvqlqp6oqpOVdU9Z1nztqp6vKoeq6qPr3ZMAFhfB/ZaUFVXJbkvyc8kOZ3kkao63t2Pb1tzKMmvJ/mJ7v5qVX3/fg0MAOtmyZn1zUlOdfeT3f1ckvuT3L5jzTuT3NfdX02S7n5mtWMCwPpaEutrkzy17fj01n3bvSbJa6rqU1X1cFXdstsTVdXRqjpZVSfPnDlzYRMDwJpZEuva5b7ecXwgyaEkb0pyZ5I/qapXvuCTuo9190Z3bxw8ePB8ZwWAtbQk1qeTXL/t+LokT++y5q+7+5vd/cUkT2Qz3gDAi7Qk1o8kOVRVN1XV1UnuSHJ8x5q/SvLTSVJV12TzZfEnVzkoAKyrPWPd3c8nuTvJg0m+kOSB7n6squ6tqiNbyx5M8pWqejzJQ0ne291f2a+hAWCdVPfOXz9fHBsbG33y5MlL8rUB4GKrqs9298aFfK4rmAHAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAwnFgDwHBiDQDDiTUADCfWADCcWAPAcGINAMOJNQAMJ9YAMJxYA8BwYg0Aw4k1AAwn1gAw3KJYV9UtVfVEVZ2qqnvOse6tVdVVtbG6EQFgve0Z66q6Ksl9SW5NcjjJnVV1eJd1L0/yK0k+s+ohAWCdLTmzvjnJqe5+srufS3J/ktt3WfdbST6Y5OsrnA8A1t6SWF+b5Kltx6e37vu2qnp9kuu7+2/O9URVdbSqTlbVyTNnzpz3sACwjpbEuna5r7/9YNVLkvxekvfs9UTdfay7N7p74+DBg8unBIA1tiTWp5Ncv+34uiRPbzt+eZLXJfnHqvpSkjcmOe5NZgCwGkti/UiSQ1V1U1VdneSOJMe/9WB3P9vd13T3jd19Y5KHkxzp7pP7MjEArJk9Y93dzye5O8mDSb6Q5IHufqyq7q2qI/s9IACsuwNLFnX3iSQndtz3vrOsfdOLHwsA+BZXMAOA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFguEWxrqpbquqJqjpVVffs8vivVdXjVfVoVf19Vb169aMCwHraM9ZVdVWS+5LcmuRwkjur6vCOZZ9LstHdP5LkL5N8cNWDAsC6WnJmfXOSU939ZHc/l+T+JLdvX9DdD3X317YOH05y3WrHBID1tSTW1yZ5atvx6a37zuauJH/7YoYCAP7PgQVrapf7eteFVW9PspHkp87y+NEkR5PkhhtuWDgiAKy3JWfWp5Ncv+34uiRP71xUVW9J8htJjnT3N3Z7ou4+1t0b3b1x8ODBC5kXANbOklg/kuRQVd1UVVcnuSPJ8e0Lqur1Sf4om6F+ZvVjAsD62jPW3f18kruTPJjkC0ke6O7Hqureqjqytex3knx3kr+oqn+uquNneToA4Dwt+Z11uvtEkhM77nvftttvWfFcAMAWVzADgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCGE2sAGE6sAWA4sQaA4cQaAIYTawAYTqwBYLhFsa6qW6rqiao6VVX37PL4d1bVn289/pmqunHVgwLAutoz1lV1VZL7ktya5HCSO6vq8I5ldyX5anf/YJLfS/Lbqx4UANbVkjPrm5Oc6u4nu/u5JPcnuX3HmtuT/NnW7b9M8uaqqtWNCQDra0msr03y1Lbj01v37bqmu59P8myS71vFgACw7g4sWLPbGXJfwJpU1dEkR7cOv1FVn1/w9blw1yT5j0s9xBqwz/vPHu8/e7z/fuhCP3FJrE8nuX7b8XVJnj7LmtNVdSDJK5L8584n6u5jSY4lSVWd7O6NCxmaZezxxWGf95893n/2eP9V1ckL/dwlL4M/kuRQVd1UVVcnuSPJ8R1rjif5ha3bb03yD939gjNrAOD87Xlm3d3PV9XdSR5MclWSD3f3Y1V1b5KT3X08yZ8m+VhVncrmGfUd+zk0AKyTJS+Dp7tPJDmx4773bbv99SQ/d55f+9h5ruf82eOLwz7vP3u8/+zx/rvgPS6vVgPAbC43CgDD7XusXap0/y3Y41+rqser6tGq+vuqevWlmPNyttceb1v31qrqqvKu2guwZJ+r6m1b38+PVdXHL/aMl7sFPy9uqKqHqupzWz8zbrsUc17OqurDVfXM2f48uTb9/ta/waNV9YY9n7S79+0jm29I+9ckP5Dk6iT/kuTwjjW/lORDW7fvSPLn+znTlfaxcI9/Osl3bd1+lz1e/R5vrXt5kk8meTjJxqWe+3L7WPi9fCjJ55J879bx91/quS+nj4V7fCzJu7ZuH07ypUs99+X2keQnk7whyefP8vhtSf42m9coeWOSz+z1nPt9Zu1Spftvzz3u7oe6+2tbhw9n82/lWW7J93GS/FaSDyb5+sUc7gqyZJ/fmeS+7v5qknT3Mxd5xsvdkj3uJN+zdfsVeeF1NdhDd38yu1xrZJvbk3y0Nz2c5JVV9apzPed+x9qlSvffkj3e7q5s/hcdy+25x1X1+iTXd/ffXMzBrjBLvpdfk+Q1VfWpqnq4qm65aNNdGZbs8QeSvL2qTmfzr4DefXFGWyvn+3N72Z9uvQgru1QpZ7V4/6rq7Uk2kvzUvk505TnnHlfVS7L5f5t7x8Ua6Aq15Hv5QDZfCn9TNl8h+qeqel13/9c+z3alWLLHdyb5SHf/blX9eDavofG67v6f/R9vbZx39/b7zPp8LlWac12qlLNassepqrck+Y0kR7r7GxdptivFXnv88iSvS/KPVfWlbP4O6rg3mZ23pT8v/rq7v9ndX0zyRDbjzTJL9viuJA8kSXd/OslLs3ndcFZn0c/t7fY71i5Vuv/23OOtl2j/KJuh9ju+83fOPe7uZ7v7mu6+sbtvzOb7Ao509wVfB3hNLfl58VfZfMNkquqabL4s/uRFnfLytmSPv5zkzUlSVT+czVifuahTXvmOJ/n5rXeFvzHJs9397+f6hH19GbxdqnTfLdzj30ny3Un+Yuu9e1/u7iOXbOjLzMI95kVauM8PJvnZqno8yX8neW93f+XSTX15WbjH70nyx1X1q9l8afYdTqDOT1V9Ipu/qrlm63f/70/yHUnS3R/K5nsBbktyKsnXkvzins/p3wAAZnMFMwAYTqwBYDixBoDhxBoAhhNrABhOrAFgOLEGgOHEGgCG+1807r3NbE+EIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xn = np.arange(data.X.min(), data.X.max(), 2)\n",
    "get_y = lambda a, b: a + b * xn\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "for i in range(models.shape[0]):\n",
    "    yn = get_y(models.a[i], models.b[i])\n",
    "    ax.plot(xn, yn, linestyle='dotted', color='grey')\n",
    "\n",
    "yn = get_y(ols['a'], ols['b'])\n",
    "\n",
    "ax.plot(xn, yn, color='red', label='OLS')\n",
    "\n",
    "ax.scatter(data.X, data.Y, alpha=.2)\n",
    "legend = ax.legend()\n",
    "ax.set_xlabel('X', fontsize=16)\n",
    "ax.set_ylabel('Y', fontsize=16)\n",
    "ax.set_title('Quantile regression with 0.05 and 0.95 quantiles')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
